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Pi erce’ 


Abstract 


s adiabatic  normal-mode  theory  of  sound  propagation  ir 


a range  dependent  ocean  duct  is  extended  to  the  case  of  a sound 
channel  with  less  than  gradual  range  dependence.  A computational 
method  is  devised  to  solve  the  ensuing  coupled  range  equations 
by  diagonalization . This  is  applied  to  a channel  with  arbitrary' 
(numerically  given)  range  dependence,  by  dividing  it  into  ranee 
segments  with  constant  properties  each.  The  local  depth  func- 
tions are  taken  as  the  A.iry  function  solutions  in  a piece-wise 
linearized  stratified  medium  with  arbitrary  sound  velocity  pro- 
file. The  method  is  applied  to  some  simple  as  well  as  to 
realistic  cases  of  range  dependent  sound  channels,  and  the  de- 
pendence of  the  mode  coupling  effects  on  the  degree  of  range 
dependence  is  determined. 
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Introduction 

Normal  mode  theory  of  underwater  sound  propagation  can,  in  its 

original  form*,  be  applied  only  to  situations  where  both  the  ocean 

floor  and  the  sound  velocity  profile  do  not  depend  on  the  horizontal 

coordinate  (range  coordinate) , so  that  the  wave  equation  may  be 

separated  into  modal  range  functions  and  depth  functions.  In  a 

realistic  ocean,  this  is  very  rarely  the  case,  but  an  approximate 

2 

method  has  been  developed  by  Pierce  in  which  gradually  range 

dependent  environments  can  be  handled  to  first  approximation.  This 

"adiabatic"  approximation,  based  on  the  Born-Oppenheimer  method  of 

molecular  physics,  utilizes  a modal  expansion  in  which  the  depth 

functions  are  taken  as  the  "local"  modes  of  a locally  stratified 

ocean  which  depend  parametrically  on  range.  Insertion  into  the 

wave  equation  then  leads  to  a set  of  coupled  equations  for  the 

range  functions,  but  these  become  uncoupled  if  in  the  adiabatic 

approximation,  the  coupling  terms  are  neglected.  The  adiabatic 

method  has  been  applied  in  our  previous  work  to  a wedge-shaped 

3 

continental  shelf  with  constant  sound  velocity  , to  a sound  channel 

4 

with  parabolic  profile  that  opens  up  linearly  in  range  , and  to 
sound  channels  with  arbitrary  (but  gradual)  dependence  of  the  ocean 
floor  and  of  the  sound  speed  profile  (which  may  depend  arbitrarily 
on  depth  also)  on  the  range^. 


No  energy  transfer  between  modes  takes  place  in  the  adiabatic 
approximation.  If  the  mode  coupling  terms  are  not  neglected,  the 
range  equations  are  coupled  and  energy  flows  between  the  modes. 


This  has  been  discussed  on  the  basis  of  coupled  power  equations 
by  McDaniel®,  and  was  considered  by  her  in  a calculation  of  sta- 
tistical  mode  conversion  effects  as  caused  by  a rough  ocean  floor  , 

For  the  deterministic  case,  the  mode  coupling  problem  has  been  solved 

Q 

in  our  earlier  work  on  the  sound  channel  with  a range  dependent  par- 
abolic velocity  profile,  by  using  diagonalization  techniques  for  the 
decoupling  of  the  range  equations. 

In  the  present  work,  similar  techniques  are  used  for  the  more 
general  case  where  the  profile  depends  arbitrarily  on  both  depth  and 
range,  and  the  ocean  floor  on  range.  This  is  achieved  by  dividing  the 
ocean  into  horizontal  layers  as  well  as  range  segments.  Linearization 
is  employed  in  the  depth  layers,  and  has  been  applied  also  to  the 

5 

range  dependent  quantities  in  the  range  segments  in  our  earlier  work 
on  the  adiabatic  case.  The  present  paper,  however,  employs  the 
approximation  of  range  segment-wise  constancy  of  the  depth  functions 
and  their  eigenvalues  as  well  as  of  the  coupling  terms,  in  order  not 
to  encumber  the  more  involved  decoupling  problem  even  further,  A 

different  approach  to  the  mode  coupling  problem  based  on  the  "slab 

9 10 

method"  , is  due  to  Kanabis  , 

Pierce's  method  will  be  implemented  for  the  coupled-mode  prob- 
lem of  an  arbitrarily  range  dependent  sound  channel  in  the  following, 
and  applied  to  simple  as  well  as  to  realistic  cases  in  order  to 
determine  the  dependence  of  the  mode  coupling  effects  on  the  degree 
of  range  dependence. 
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I.  Local  Eigenmodes  of  the  Wave  Equation 
We  seek  a solution  for  the  time  harmonic  pressure  field 

-=> 

p (r)  of  a point  source,  located  at  the  position  rQ  in  an  inhomo- 
geneous ocean  of  density  (r)  and  propagation  constant 
k K ) - UJ  /c  ( -f  ) , as  determined  by  the  wave  equation 

$(■?)*  C'r)V,p(r)^-i-WL  ( r*  ) - J fr-  r0  ) . (la) 

Adopting  a cylindrical  coordinate  system  with  the  z-axis  pointing 
vertically  downward,  we  shall  assume  the  sound  speed  c(r)  to 
depend  arbitrarily  on  both  the  depth  z,  and  on  the  horizontal 
range  coordinate  ■&  (x,  y)  . If  the  medium  (which  includes  the 
ocean  bottom)  is  layered,  the  boundaries  of  the  / th  layer 
z « z ~ may  a^-so  depend  on  ^ while  the  ocean  surface  z = z^ 

is  flat.  Densities  <p  (z)  are  taken  constant  within  each  layer; 
they  v/ere  chosen  equal  to  <p^  in  each  water  layer  but  could  be 

set  equal  to  a different  constant  in  each  bottom  layer.  Defining 

■ — * 

a velocity  potential  <p  (r)  such  that 

^ (*?  ) = - V Cp  (?)) 

jj  r-P)  = P (?)  ?$(?)/&  (lb) 


(v  being  the  particle  velocity) , one  has  in  each  layer  a 
Helmholtz  equation 

(?)  + kx{*,t  )$(■?)*  ?(■?-£).  (lc) 

2 

Following  Fierce  , we  attempt  a solution  in  the  form  of  a quasi- 
separated  normal-mode  sum 
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'n 


(2) 


in  which  the  "local  depth  functions'1  -U,^  (7^  <p  ) satisfy  at  each 
range  point  ^ (considered  as  a parameter)  the  depth  equations 
pertaining  to  a stratified  medium. 


9 (-£,9 ) 

gl1 


K! 


U,?)  Uyv(*/?>  = 0. 


(3) 


The  eigenvalues  k.^  ) of  the  local  modes  enter  through  the  local 

vertical  wave  number 


Mx~ 


(4) 


of  the  nth  mode.  They  are  determined  by  the  boundary  conditions 
satisfied  by  the  local  depth  functions  as  used  in  Eq.  (2) , 


fj.  ^Ltv  / f ) ~ j ^ ,?)  , 

[9m  ./Qz]  = Q 2^ /3z  ] 


(5) 


at  the  boundary  between  layers  Ji  and  £+1  ( ((?  ) being  the 

boundary  as  approached  from  the  ,/th  layer,  and  (%  ) as 

approached  from  the./f+lst  layer) , corresponding  to  continuity  of 
pressure  and  normal  particle  velocity.  Inserting  Eq. (2)  into 
Eq.  (lc)  and  using  the  orthogonality  of  the  depth  functions  which 
follows  from  Eqs.  (5), 


r r,p 
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/ Lf  U,p  ) ch  ^ £ 


n »v  > 


(6) 


we  obtain,  following  our  earlier  derivation^,  the  coupled  system 
of  (partial  differential)  range  equations 


(7) 


in  which  the  coupling  coefficients 


M (i 

yi  *yw  I 


M 


CO 

) = / [ f fe)/fw']  P ) Vf  Uhv  (55,(3  ) di 


(8a) 


(8b) 


couple  the  nth  and  the  mth  inodes.  We  here  used  the  notation 

V0  = (d/dx,  2/2y,0),  and  took  the  source  position  as  r = (z  , 
y o o 

« 0 ).  The  derivation  of  Eq.  (7),  due  to  our  use  of  Eqs.(5), 

holds  also  for  the  present  case  where  sloping  layer  boundaries  are 
admitted.  However,  Eqs.  (5)  will  guarantee  the  boundary  conditions 
of  the  total  acoustic  field, 

[ty ,n  -[Vf  /9n.]  t 


(9) 
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to  be  approximately  satisfied  only  for  gradually  sloping  boundaries 
and/or  small  relative  density  differences  between  layers}^  i.e., 

r (fui  - « i no) 

(where  CX  is  the  angle  between  the  boundary  normal  and  the  z-axis) , 

as  can  easily  be  shown.  This  will  be  assumed  in  the  following. 

Note,  however,  that  no  restrictions  are  imposed  on  the  range  gra — 

dient  of  the  velocity  profile. 

The  coupling  terms  of  Eqs.  (8)  may  be  recast  in  a simpler 
3 

form  which  shows  their  origin  from  the  range  dependence  of  the 
sound  velocity  profile  and  of  the  boundaries,  res  pectively,  by 
employing  partial  intergration  and  using  Eq.  (3),  with  the  following 
results.  For  the  coupling  terms  which  are  small  of  first  order 
in  the  range  dependence,  one  has 
— ^ 

"t  'C.r  (?)  (■»■>■  w)  (lla) 

(p)  = (?)  (11b) 

with  the  volume  contribution 

r?)/fX  (?)  - ^ (?)],  (12a) 
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— * 

^'KW.  (?)  = J If  (12b^ 


and  the  surface  contributions 


O’  (f)  = n„‘”  (f)/[l<:  ■p’J-Cif )], 


(13a) 


^-u»v  ^ = il  - lL  - ?)  (13b) 

J~J  L <P  1 v 


*iV) 


where  L is  the  number  of  layers,  and 


(13c) 


further. 


<■ (s) 
oarx 


<p ) 


fft/pw  )[<(*,*)%  i\  „ . 

x-i  ■ f ( p J 


(13d) 


For  the  second-order  coupling  coefficients,  we  shall  only  quote 
the  results  for  the  volume  contributions: 


O ' P > = N‘i 


), 


(14a) 


NZ(?>  = - jz  o <p>-Wf>- 


(14b) 


f>'+  V-j  rvu 


The  surface  contributions  are  absent  if  the  boundaries  are  flat. 


r-'T 


II . Solution  of  the  Coupled  Range  Equations 

In  our  previous  work  on  mode  propagation  in  a range 
„ , . 3-5,8 

dependent  environment  , two  important  special  cases  were 

considered  in  which  range  dependence  occurs  either  in  the 

Cartesian  x-direction  only  (x-case) , or  in  the  ^-direction  of 

cylindrical  coordinates  only  ( ^ -case) . For  sufficient  distances 

from  the  source,  the  results  of  these  two  assumptions  were  shown 
4 

to  coincide  . We  shall  here  restrict  ourselves  to  the  ^ -case 
only,  so  that  } and  introduce  the  new  range 

function  f ( <o  ) via 

n S 


= f'1'*  f.  (?) 


(15) 


Inserting  in  Eq.  (7)  leads  to  the  coupled 
system  of  ordinary  differential  equations 

I*," - (**> 

= 17 ^ t^.=) /<pv^  1 f'2xp‘/'!)  J/pj 


where  M'  ( <p  ) is  defined  similarly  to  M’  in  Eqs.  (11) -(13) 
but  with  replaced  by  d/d<=  , and  where  the  further  coupling 

terms  are 


(?)  “ M <?)  - f?)- 


(16b) 


4 -2 

Note  that  as  shown  before  , the  term  (2  o ) f in  Eq.  (16a)  is 

j n 

negligible  for  distances  from  the  source  in  excess  of  ~1  km,  for 


frequencies  up  to  several  Hz.  Similarly,  the  second  term  in  Eq.  (16b) 
is  negligible  as  compared  with  the  second-last  term  in  Eq.  (16a) 
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in  an  even  wider  region.  We  write  Eq.  (16a)  in  matrix  notation  as 


(16c) 


where  f(£  ) is  a column  vector  with  elements  f^^j) « and  k (y)  a 


/ —r 

diagonal  matrix  with  diagonal  elements  k^  (p) ; U is  the  column 


vector  of  the  source,  with  components 


= [f  (^o)/pvv]  Uk^)  ^ J'/f)  .u^  (Z-o;0).  (16d) 


Our  goal  will  be  to  devise  a transformation  which  uncouples 


the  system  of  Eq.  (16c),  at  least  in  an  approximate  manner.  This 


will  be  accomplished  by  first  dividing  our  range  ^ into  a finite 


number  of  inter  als  with  constant  vertical  boundaries,  across 


which  <^>(z )<p  and  df/dQ  have  to  be  continous.  In  the  approximation 


of  Eq.  (10),  this  is  shown  by  integrating  over  z and  using  ortho- 


normality of  u^,  to  lead  to  continuity  of  f^  and  dfn/dp,  i.e.  to 


(pT  ) = ffiT, ) 


V‘  ?<*: 


where  is  the  outer  boundary  of  the  t’-(&  range  interval, 

coinciding  with  the  inner  boundary  C,  of  the  it i s<  range 

1 -t  1 1 


-n 


J 


interval . 
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In  our  earlier  study  of  the  "adiabatic"solution  of  Eq.  (37) 
where  the  coupling  terms  were  neglected,  we  approximated  the  quan- 
tity ICl?)  ♦!/*»?*  in  each  range  by  the  linear  expressions 
d^rT)  , joined  continuously  at  each  segment  boundary  where  the 
actual  values  were  adopted,  except  in  the  first  (source)  interval 

and  in  the  last  one  (reaching  to  infinity),  where  the  constant 

2 2 2 
values  k^  (o)  and  k^  (<»  ) , respectively,  were  adopted  for  k^  (^>) . 

This  leads  to  Airy  function  solutions  in  the  general  intervals, 

and  to  Hankel  functions  H in  the  two  extreme  intervals,  which 

2 

become  trigonometric  solutions  if  the  term  l/4^>  is  neglected. 
This  procedure  would  also  be  possible  in  the  present  case;  the 
coupling  terms  would  then  have  the  form,  e.g.. 


(18) 


2 

with  depending  on  p only  through  u and  u , Since  k (z,p) 

'"Tit  ) TV  m ' ) 

is  also  being  approximated  by  linear  expressions  C(z)p4-D(z)  in 
order  to  represent  the  experimental  range  (and  depth)  dependent 
sound  speed  profile  which  is  given  to  us  only  at  a finite  number 
of  range  points  (chosen  as  our  range  interval  boundaries).  We 
shall,  however,  for  reasons  of  greater  simplicity  approximate 
in  each  range  interval  the  eigenvalues  k^^p)  and  depth  functions 
<4.  (z.?).  and  hence  also  the  coupling  terms,  by  their  values  at 
the  midpoints  of  each  range  interval,  insuring  accuracy  by 
introducing  a finer  subdivision  into  range  intervals  if  necessary. 
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We  are  then  able  to  uncouple  Eq.  (16c)  by  matrix  transformations 
involving  constant  transformation  matrices  in  each  range  interval. 

g 

Following  our  previous  approach  used  in  the  treatment  of  mode 
coupling  for  a range  dependent  parabolic  profile,  we  first  set 

(?)  = ^ ^ (f)  (19) 

which,  with  a constant  matrix  leads  to 

[dV<v+  K^jf -5/  u-ixr'-wr 


where  we  neglected  l/4(?  , and  called 

Kz(?)  = X k SV  (21a} 

2 = J?"1  M'  Sx  1211 

W = 1 V ^ ■ (2it 

The  matrix  will  be  so  chosen  that  X is  diagonal,  with  dia- 


(21b) 


(21c) 


gonal  matrix  elements  that  will  be  called 


X . (Tr 


(This  is  possible 


with  a constant  matrix  3 since  as  mentioned  above,  M’  is  taken 
constant  in  each  range  segment) • Next,  the  derivative  terms 
r which  now  multiply  the  diagonal  matrix  X.  t may  be  elim- 
inated by  standard  methods.  We  make  the  transformation 


r - 


(22a) 


where  (Q  is  taken  to  be  a diagonal  matrix  with  elements  CT  . 

'n. 
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Substituting  Eq.  (22a)  into  Eq.  (20)  and  equating  the  coefficients 
/ 

of  to  2ero#  we  find 

; (22b) 

'TV  ) 

constant  factors  in  Eq.  (22b)  may  all  be  chosen  equal  to  unity. 

The  transformed  range  equations  are  still  coupled  (although 
without  first  derivatives) , and  have  the  form 


-A2/  + (K2’-*-'"/  )y  - 1 U . 


(23) 


The  source  term  is  unaffected  by  the  transformation  of  Eq.  (22a) 
since  U contains  a delta  function  in  (p  , and  = is 


the 

identity 

matrix. 

The  elements  of 

the  non-diagonal 

matrices 

"2 

r** 

K 

and  W are 

given  by 

KL 

^ (A  yt  - )<p 

K 2’ 

' 'm  >rt.  I 

(24a) 

w 

SY\  >w 

^ ^ M ) (p 

In  the  spirit  of  our  approximation  procedure,  these  will  again 
be  replaced  by  their  values  at  the  midpoints  of  the  p-intervals, 
rendering  the  expression 

T = K 1 + V/  - A1  (25) 
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where  is  chosen  such  that 
2 


- i 


T 51  s A 


(28) 

is  a diagonal  matrix  with  elements  * This  has  finally  brought 

us  to  the  uncoupled  system  of  linear  equations  with  constant  coeffi- 
cients > — a. 

d.zP/oL9"+  AP  = U * , 


where  the  new  source  term  is 

— * . 

^ \ “ -L 


U = (S.SJ--  u . 

Its  general  solution  is  in  the  source-free  region: 

N 

T ~ H (cx^  -h  pt- 

- 1 ^ 


(29 


(29b) 


(30a) 


where  N is  the  total  number  of  modes. 


^ a 1/z 


(30b) 


7ZT ^ are  the  basis  vectors  with  components 

- <*>=> 

and  are  2N  arbitrary  coefficients.  If  we  call  the  overall 

transformation  matrix 

U = 3^  , (31a) 

our  original  range  functions  of  Eq.  (15)  are  obtained  as 


{ = U?  , 


(31b) 
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with  components 


Itv 


(31c) 


The  Green/s  function  of  Eq.  (29a)  which  describes  the  behavior  of 
the  solution  at  the  source  can  be  obtained  from  Eq.  (31c)  in  the 
well-known  way,3  5 see  below. 

As  the  last  step,  we  shall  have  to  carry  out  a matching  of 
the  range  functions  and  their  derivatives  at  the  range  interval 
boundaries,  see  Eq.  (17).  In  the  last  interval  which  reaches  out 
to  infinity,  we  assume  no  range  dependence  so  that  k^ (p) =constant 
=kp(  oo  ),  and  have  the  solution5 


v1'1  H“* 


P V ' r • ' (32a) 

(no  mode  coupling  takes  place  here) , comprising  only  outgoing  waves 
with  undetermined  modal  amplitudes  A^.  In  the  general  (ith)  interval, 
we  adopt  Eq.  (31c)  with  superscripts  i on  f ) 0(  t (5^  an<^ 

w ; in  matrix  notation,  this  gives 


r?)  = Ui  cT‘  C<?) 


(32b) 


where 


is  a column  vector  with  components 


(32c) 
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Introducing  the  diagonal  matrix  A with  diagonal  elements  Ap , 

- ^ 

and  the  vector  H(^>)  with  components  (o")=  <p1/2-  H,,'^ 

Eq.  (32a)  reads  in  matrix  form 

?).  (32d) 

If  we  designate  the  last  range  interval  before  the  infinite  interval 
by  i=M,  the  matching  conditions  are  (at  the  boundary  9=f>M): 


(J|V|  = A H (pM) 

UM  ~ AH'  ip*)- 

We  introduce  the  two-component  vectors  (2D  space) 


Ct. 


M 


CK 


M 


1 


H 


rt 


H 


(fn>  = 


H),  tf„) 

Hp  <f«)/ 


and  the  2x2  matrix 


if  = 


9%  ?n 


ii  Hi  ws^PAi. 


with  inverse 

Then,  Eqs.  (33)  are  solved  by 


(33) 


(34a) 


(34b) 


(34c) 


(35a) 


The  notation  may  be  further  compacted  by  intooducing  hypervectors 
**  M 


oC  j f-j  (9n)  whose  components  (each  corresponding  to  a mode) 
are  the  2-component  vectors  c*  ^ ^ f j (?m)  • respectively,  as  well 
as  diagonal  hypermatrices  f whose  diagonal  elements  are  the 


2x2  matrices  ifL  (so  that  the  elements  of  the  diagonal  hyper- 
~ l * 

matrix  '4  are  the  2x2  matrices  ) . Eq.  (35a)  then  reads 

o(M  = ('if  ^ Un*  A H ) , (35b) 


where  the  previous  mode-space  matrices  U ^ and  A are  understood 

as  being  diagonal  in  2D  space. 

Applying  Eq.  (17)  for  the  match  between  intervals  i and  i+1, 
we  have 


ir  r = u4rl  } ip<) 


(36) 


which,  in  our  hyperspace,  is  solved  by 


)_i  u/  1 U**1  if 


<£-> 

cx 


(f 


(37) 


Finally,  we  match  at  the  boundary  between  invervals  i=l  (source 
interval)  and  i=2.  In  the  former,  there  is  no  coupling,  and  the 
solution  is"* 

f,  <?' = a;  n:\Mpb  p;  /Ok,™?)] 


(38) 
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where  A ° are/known  amplitude  factors  which  we  again  arrange  in 
a diagonal  matrix  A ^°\  Introducing 


U (i,l>  , X i/A  1 1 x 

H p ) s 9 Wo  ( 


(39a) 


and  forming  the  known  matrix 


^ - 


H»(l1  (px) 


H,'11'  (fi) 


\ 


we  obtain  from  Eqs.  (17): 


(39b) 


A"”  t'1  Ui r "'Pi-  oi3- , 


where  hypermatrix  notation  was  introduced  as  before.  In  this  way, 
all  coefficients  O^'1  and  ^3^*  have  been  determined;  they  are 
in  the  ith  interval: 

2c  * - u/4  IT*1  K.*)'1  u.^uifZ 

r . . (£" ft. )' ■ 1 UM- 1 A H (Pn ) , («-) 


and  in  the  source  interval: 


A1’1  S 1 = "fr1-  i'e2'f^)~L  ux'‘  ui  e 


( ‘.ft.-. V/  IT  (•€' O H (P„) 


(41b) 


These  expressions  still  contain  the  unknown  coefficients  Av  and  A, 
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but  as  was  shown  earlier,  matching  to  the  source  leads  to 


(o) 


Alp 


HHo/pl  - p>pL ) f 


(42) 


w 


so  that  in  principle,  the  remaining  unknowns  A p (the  outgoing 
mode  amplitudes  at  infinity)  can  be  determined  by  the  require- 
ment that  Eqs.  (41b)  and  (42)  be  satisfied.  The  implementation 
of  this  prescription,  however,  requires  an  iterative  approach 
which  will  be  described  below. 

III.  Numerical  Results  and  Discussion 
The  above  procedure  will  now  be  applied  to  the  calculation 
of  transmission  loss  in  some  examples  of  channels  with  simple 
and  also  with  realistic  velocity  profiles.  In  order  to  be  able 

to  handle  the  la.tter  case,  we  shall  use  the  previously  employed 
5 

method  for  calculating  local  depth  functions  tt^fz,^)  and  their 
2 

eigenvalues  k^  (q)  in  a numerically  given  profile,  which  consists 

in  dividing  the  ocean  into  (horizontally  bounded)  depth  segments  or 

2 

layers,  in  each  of  which  k (z,  (>  ) is  approximated  linearly,  so  that 

the  depth  functions  are  given  by  Airy  functions.  The  layer  below 

the  (locally  flat)  ocean  floor  is  simply  modeled  as  a liquid 

with  constant  bottom  density  , and  with  a H-independent  local 

sound  velocity  /Cg(^)  , which  reaches  down  to  z . The  eigen- 

2 

values  k^  (<£)  are  then  obtained  by  satisfying  the  boundary  conditions 
of  a pressure  release  surface  at  the  ocean  surface  z=zgr  i.e.. 


; i 
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The  following  computer  program  is  based  on  the  foregoing 
theory;  it  computes  range  functions  and  solves  the  coupled- 
mode program  of  sound  propagation  in  a range  dependent  ocean. 

The  main  program  GPC  (pp.  24  - 32)  computes  the  propa- 
gation loss  as  a function  of  range  and  depth.  It  calls  all  the 
subroutines;  it  calculates  the  range  functions,  it  calculates 
the  coefficients  of  the  linear  combinations  of  the  range  func- 
tions in  all  the  range  segments. 

The  subroutine  NAGL  (p.33)  interfaces  the  range  function 
and  depth  function  programs;  it  calls  the  depth  function  program 
which  is  described  in  Volume  I of  this  report. 

The  subroutine  FORM  (pp. 34-37)  calculates  the  U1  s and 
the  q-eigenvalues 

The  subroutine  C0EF1  (pp. 38-40)  calculates  the  mode  coef- 
ficients in  the  source  segment  and  in  the  infinite  segment  by 
matching  directly  (since  all  the  internal  boundary  conditions  are 
here  lumped  together  into  one  matrix) 

The  subroutine  HMDIAG  (p.41)  diagonalizes  an  antihermitean 
matrix  (by  relating  it  to  the  diagonalization  of  a Hermitean  one) 

The  subroutine  MXEL  (p.42)  computes  integrals  J unumz  dz 
containing  Airy  functions 

The  subroutine  HANK  (p.43)  calls  the  Hankel  function 


subroutine 


The  subroutine  HANGL  (p.43)  computes  the  Hankel  functions 


of  argument  X 

The  subroutine  DAIRY  (p.44-50)  is  the  same  as  in  the 
depth  function  program  (Vol.  I of  this  report) 

The  subroutine  PLT  (p.51-52)  is  a plot  routine  to  plot 
the  propagation  loss 

The  last  page  (p.53)  tells  what  has  to  be  loaded  into 
the  range,  depth,  and  plotting  programs. 
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ty  ar  c 


COMMON  /DEPTH/  XKMII5.4)/ 
1U0C15)  1 ' 


COMMON  /FORMC/  XT'  ^<4)  » ZB ( 4 ) ,-NM » NL  , NLM , NLP » XK  ( 4 » 4 ) , PEN  < 4 > 
1»XNAG(4) 

DIMENSION  XKT(-’;,  .<i5'4>’ 

lU<15»15)»UI(l.J,i5),U<l5> » AM  (30  >30)  , V(  IS)  .BUM (30 ,30)  , 

2A  ( 15) » AO  ( IS)  (3,  15)  »BET  ( 3 , IS ) » US (2, IS.  IS)  . UIS(2»15.15)  t 

30  S ( 2,  15 ) »Fi\  ( 15  ) r j,V  \ 15  ) » IOUT ( 15 ) » G < 30 ) , XKNT  ( 15 ) » N.MOE'E  \ 4 > 
COMPLEX  L'»  01  . Q / 1'.'1 » SUM » E YE  > V . DUM » A , AO  , ALF  . BET  » US » UIS » QS » SUM2 » 
1FR » UN  »Ci»wi,m<HlP»  H2 . H2F-  > ill » A2 » A3 » A4 , A5 . TEM  . OR 
DIMENSIONED*  15>  » UF<D  < 1 5 » 4 ) * XKND  < 1 5 » 4 > 

DOUBLE  PKLClSlON  XNAG » UOD » URD » XKND 
EYE=<G.,i7> 

PI=3. 1415Y 
TYPE  720 


!15.15),UI(1  . 15)  » <<  i5  > ’ AIK  30 » 30) » V( 13) > BUM (30.30) , 

; 15), AO  US), ■ (3, 15) » DET (3.15)*US(2,1S. IS) ,UIS(2.15,15 
> < 2, 15 ) ,FR  < ;5  J.j'15)  » I OUT « 15)  r G(  30)  , XKNT  (15),  NMOLE  ( 4 > 


-V  ' ‘ 


-A  ; 

720 

FORMAT!'  TYPE  0 FOR  NO  COUPLING,  1 FOR  COUPLING') 

. © 

i 

i 

ACCEPT  2 , I COUP 

TYPE  1 

• ; 

j 

1 

FORMAT ('  INPUT  NUMBER  OF  MODES,  LAYERS,  AND  SEGMENTS') 

ACCEPT  2»NM»NL»NS 

. • 

i 

TYPE  2, nm.nl, NS 

• 1 

I * 

n 

FORMAT <5G) 

NLM=NL-1 

• 

! 

NLP=NL+1 

• 

NSM=NS-1 

• 

i 

NSP=N3+1 

TYPE  460 

i 

: e- 

• . 

! 

460 

FORMAT  ('  INPUT  NUMBER  OF  MODES  FOR  EACH  F’RCFILE ' ) 

ACCEPT  2, (NMODE(I) »I=1,NSP) 

TYPE  2 , ( NMODE ( I ) , J=1 »N3P ) 

e ; 

TYPE  6 

»; 

6 

FORMAT ('  INPUT  FREQUENCY') 

ACCEPT  2,F 

• 

^ • 

TYPE  2 »F 

© 

F=2.*PI*F  ' 

• 

TYPE  700 

700 

FORMAT < ' INPUT  CONVERGENCE  CRITERION') 

ACCEPT  2.FRAC 

TYPE  2,FRAC  _ 

■ • 

9 , 

TYPE  3 

•© 

fi  0 

H 

■1  e 

«< 


, 


3 FORMAT ( ' INPUT  SEGMENT  POSITIONS') 

X ( 1 ) =2000  < 

X (2  >=20000 . 

X ( 3 > =40000 • 

TYPE  2t <X(I> ,I=1»NS> 

ACCEPT  2»(X(I)»I=1* NS ) 

XNAG ( 1 ) =0 . 

DO  200  1=1 » NS 
200.  XNAG(IH) --X(I) 

TYPE  4 

A FORMAT ( ' INPUT  SOURCE  DEPTH.  RECEIVER  DEPTH,  NL-1.  PROFILE' , 
1'  DEPTHS,  NS f 1 BOTTOM  DEPTHS') 

ACCEPT  2 » ZS » ZR . < Z ( I > » I =2  » NL ) » < ZB  < I > ,X*1,HSP> 

TYPE  2,ZS,ZR , ( Z < I ) , I=2,NL  ) , (ZD  ( I ) > I -1 » NOP-) 

Z<1)=0. 

DO  101  IZS=1,NL 
IF  < ZS . LT • Z ( IZS  H ) ) GO  TO  0 
101  CONTINUE 

8 TYPE  5 • 

5 FORMAT < ' INPUT  SOUND  SPEED  PROFILE  FOR  EACH  SEGMENT') 

DO  002  IS---1  ,NSP 
802  XK ( 1 , IS) =1520 « 

XK(2, 1 )= 1 5 1 5 . 

XK  ( 2 » 2 ) = 1514 . 5 
XN( 2,3) *1510. 


I 

4 


< 


■<  o 

I 

•{ 

» 


© 


ji  °l 

i 

* © i 


9 

i # 

i 

r 


XK  (2,4 ) =1000  • 

XK(3, 1 )=1490. 

XK<3,2)=1490. 

XK(3,3)=1490. 

XK (3,4 ) = 1490 • 

XK (4*1 ) = 1540  « 

XK(4f2>=1540. 

XK ( 4 , 3 ) = 1540 . 

XK(4,4)=1540. 

DO  100  IS=1 , NSP 

C ACCEPT  2 » < Xl\(  ILfIS)fIL=1fNLP> 

TYPE  2f(XK(IL,IS)fIL=1fNLP) 

DO  100  IL=1 fNLP 
XK(ILfIS)=(F/XK(ILf IS)  )**2 
100  CONTINUE 
TYPE  7 

7 FORMAT ( ' INPUT  RELATIVE  DENSITIES  FOR  EACH  LAYER' ) 

DO  803  IL=1,NLP 
803  DEN ( IL ) = 1 . 

TYPE  2, (DEN(IL) ,IL=1,NLP) 

C ACCEPT  2, (DEN(IL)fIL=1fNLP> 

TYPE  40 

40  FORMAT ('  TYPE  NUMBER  OF  MODES  TO  BE  OUTPUT') 

ACCEPT  2.N0UT 

TYPE  2,N0UT 
N0UT2=2*N0UT 
IF(NOUT.EO.O)  GO  TO  41 
IF(NOUT.EO.NM)  GO  TO  43 
TYPE  42 

42  FORMAT ( ' TYPE  WHICH  MODES  ARE  TO  BE  OUTPUT') 

ACCEPT  2f(I0UT(I)fI=1 fNOUT ) 

GO  TO  41 

43  DO  44  1=1 fNM 

44  IOUT < I ) = I 
C 

C DETERMINE  DEPTH  FUNCTIONS 
C V 

41  CALL  NAGL ( NM  , NLP  f NSP  f NMCDE  f Z f ZB  f XK » XKN  f UO  f UR  f ZS  f XKND  f UOD  f URD ) 

C 

C FORM  PRODUCT  Oc  MATRICES  ON  PACE  44 . . .U2*H2 ( XI > *H2I  C'2 > *U2I*U3* . . . 
C FOR  USE  IN  CONNECTING  SOURCE  INTERVAL  WITH  INFINITE  INTERVAL. 

C SIZE  OF  MATRIX  IS  NOW  TWICE  THE  NUMBER  OF  MODES  TO  ALLOW  FOR 

C INGOING  AND  OUTGOING  SOLUTIONS.  MATRIX  U IS  DIAGONAL  IN  THE  2 

C SOLUTIONS.  MATRIX  SCRIPT-C  IS  DIAGONAL  IN  MODE  NUMBER. 

C 

‘ NM2=2*NM 
DO  106  I=2fNM2 
JM'-I-l 

DO  107  J=1fJM 
AM(IfJ)=0. 

107  AM(J,I)=0. 

106  AM( I » I >=1 . 

AM(1f1)=1. 

C TYPE  201f<(AM<IfJ)fJ=1fNM2)fI-1f,<,'1:’> 

201  FORMAT ( ' 201  AM'/(12G>) 

DO  108  IS=2fNS 
ISM  IS- 1 
ISTEM=IS 

CALL  FCRM<ISTEMFUFUlFaFlCCL'P?NMr-  -> 

C TYPE  621f(G<I)fI=1fNM) 

621  FORMAT  <'  621  Q'/6G) 

DO  113  1=1 fNM 
DO  ll'J  J--1  fNM 
LC(ISMfIfJ)  I'(IfJ) 
uc  u.^cieMfIf.))  un i f j) 

DO  ii 9 1=1, NM  ... 


23HIS  PACE  IS  BEST  QUALITY  PRACTICABLE 
EROU  COPY  EURKISHED  TO  L^Q  

O , 


0 


O 


*> 


o 


o 


o 

9 


o 


* 


o n n n n n n o noono  oonn  onn 


- 2€»- 


ii?  us<ismfI>=o<i> 

MULTIPLY  FROM  RIGHT  BY  U 

DO  109  1=1 fNM2 
no  109  K=1,NM 
no  109  J=1f2 
3UM-0. 

no  110  L=1fNM 

110  SUM  3U.!HAHCI,2*<L-1>+J>*U<L,K> 

109  BUM-C  I » 2>Jt < 1<—  1 >+ J)=SL'M 

TYPE  222 f ( ( BUM JI»J)*J=1 fNM2)  » 1=1 »NM2> 
>5  FORMA T( ' 222  BUM'/<12G)) 


MULTIPLY  FROM  RIGHT  BY  SCRIPT-C  EVALUATCB  AT  LEFT  OF -SEGMENT 

DO  111  L-l fNM2 
BO  111  1=1 »NM 
11=2*<I-1> 

C1=CC03 (0 ( I ) #X ( ISM) ) 

S1=CSIN<CKI)#X<I3M>> 

AM(L»Il+l)=nUM(L»Il il>*Cl-BUM<L»Il+2>*Q(I>¥Sl 

111  AM(Lf  I1+2)=BUM (L » 11  + 1 >*31+nUM(L » Il+2)*fU  I )>!<C1 

TYPE  223  f < (AH(IfJ) t J=1 »NM2) » 1=1 ,NM2) 

223  FORMAT  ('  223  AMV(12G>> 

MULTIPLY  FROM  RIGHT  BY  SCRIPT-C  INVERSE  EVALUATE!!  AT  RIGHT  OF 
SEGMENT 

BO  112  L=1*NM2 
no  112  1=1 fNM 
I1=2*<I-1) 

C1=CC03<0( I >*Xf IS) ) 

S1=CSIN(Q<I)*X<IS) ) 

BUM (L  f 1 1 + 1 ) =AM < L » 1 1 +1 ) *C1 f AM (LfI1+2)*S1 

112  BUM <Lf 11+2) =<-AM<Lf  II HO  »SHAM  (Lf  II I 2)#C1  )/CJ<  I ) 

TYPE  224 / < ( BUM < I f J)  f J=1  »N’M2)  1 1=1 1 NM2 > 

224  FORMAT ('  224  BL'M'/<12G)) 

MULTIPLY  FROM  RIGHT  BY  U INVERSE 

BO  113  1=1 fNM2 
BO  113  K=1.NM 
BO  113  J=1f2 
SUM'-O . 

BO  114  L=1 • NM 

114  SUM~SUM-KOUM  ( 1 f 2*  (L-l ) + J> *UI  (L f K ) 

113  AM  <1  / 2 1:  ( K- 1 > v J > =SUM 

TYPE  225 /<<A?KIfJ)f J"1 t MM2 >.I=1» NM2) 

225  FORMAT < ' 225  AMV(12G)> 

103  CONTINUE 


• 

% 

% 

% 

% 

I 

« 

I 

%> 

%) 


MULTIPLY  FROM  LEFT  BY  SCRIPT-II  INVERSE  * O 

X2=SCJRT (X(  1 > ) 

BO  115  1=1 fMM  O 

I1=2*<I-1> 

RHO- SORT < X\N < I , 1 > > ( 1 ) 

CALL  HANK <RHO f 1 11  • HI P f H2 f H2P ) O 

TYPE  243  f I fR'! !0  fHI  fHIP  f 1 12 f! I2P 

XNTEM=RH0/X(1> 

A1  = X2  Till  * 

A2--X2SH?  • » ... 

A3--XK  rrM-.;<x2*H  1PHU/C2.  TX2  > 

rc;i;  xr.  *;  I2P  JH2/<2.  ■'  x:* ' * 

A5=Ai#A4-A2*A3 


• w 
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THIS  PAGE  IS  BEST  QUALITY  PRACTICABLE 
roOM  COPY  FURatlSHED  TO  PDO  ^ 


TE/I-Al 

A1^A4/A5 

A4-TEM/A5 

A2=-A2/A5 

A3=-A3/A5 

TO  115  L=1 » MM2 

DUN  ( 1 1 f 1 » L ) -A1  *AM  < 1 1 + 1 , L ) f A2YAM  < 1 1 +2 » L ) 

115  DUM(Il+2»L)=A3*Ai-i(IlHfL)+A4:XAM(Il+2»L) 

TYPE  226  f ( ( HUil  < I , J > r J-l  t NM2 ) » I~1 , NM2 ) 

226  FORMAT ( ' 226  0UMVU2G)) 

MULTIPLY  FROM  RIGHT  BY  H 

XNS=SQRT  < X < NS  > > 

HO  116  J-l r NM 

ji=2*(j-n 

RHO=SORT  < XKN ( J , NSP > > *X ( NS ) 

CALL  HANK ( RI-IO > Hi  » HIP » H2 » H2P,) 

TYFE  243 » J » RMO  * HI > HIP » H2 » HUP 
XKTEM=RHO/X  < NO ) 

DO  116  I---1.NM2 

116  AM!  If  J)-DUn(I»  J1  + 1)*XNS*H1+DUM«I»  Jl+2) * (XNTEM#XNS*H1P+H1/ (2 . * 
1XNS ) ) 

TYPE  227 » < ( AM ( I » J > » J~ 1 » NM ) » I = 1 » MM2  > 

227  FORMAT ! ' 227  AM'/!6G)> 

FORM  MATRIX  TO  BE  USED  IN  DETERMINING  COEFFICIENTS  IN  SOURCE 
SEGMENT 

DO  117  I-l f NM 

117  V!I)=UO!I>*DEN(IZS>/(4.*EYE> 

TYFE  223 » < V ( I > » I = 1 * NM ) 

228  FORMAT  ( ' 229  W(6Q>) 

FIND  COEFFICIENTS  FOR  SOURCE  INTERVAL 

V 

• CALL  COEF1 <NM» AM»V»A»AOf ALFtBET rFRAC) 

DETERMINE  COEFFICIENTS  FOR  CTHER  INTERVALS  USING  RELATION  OH 
. PAGE  44 

DO  120  I=2»NM2 
JM=I-1 

DO  121  J=1 » JM 
AM ( I » J) ~0 • 

121  AM(JfI>=0. 

120  Ai!  ! I » I > *1  * 

AM ( 1 • 1 ) =1. 

TYPE  2 40 1 ( ( AM ( 1 1 J ) r J - 1 » NH2  > • I - 1 f NM2 ) 

240-  FORMAT! ' 240  AM'/( 120) ) 

MULTIPLY  FROM  RIGHT  BY  SCRIPT  C INVERSE  EVALUATED  AT  BOUNDARY  OF 
INFINITE  SEGMENT 

DO  122  L=1»NM2 
DO  122  I-l *NM 
Il=2t ( I - 1 < 

C1~CC03<C1S<NSM»  I )*X(NS) ) 

SI  COIN ! CJS  <N3M * I >1X( NS  ) ) 

DLM ! L / 1 1 ( 1 ) -AM  <L  » T 1 H > tr  1 1 1+1) tGl 

122  DUNiLf Iif2>-< -AN (L. II 4 1 > ' F : +A*(L|  H+2  ) *C1  > 'OS (NSM  • I ) 

T'.PE  241 1 < (DUM(1 1 J)  f J It.- Ml)>I  = I»NM2> 

241  F ORMAT  < * 241  DUN ' / ( 120 >> 

ii’.i-i  ;i  i.y  f rom  r.-r;  ' by  u inverse 


! C 

e 

i 

! C 

i 

• r 

r 


/ 


c 

c 


SEGMENT 


XNS-SCRT (X ( N2  > > 

Hu  125  J=  1 /NM 
( J-l ) 

RH0=S3RT < XKN < J / MSP ) ) *X ( MS  > 

CAUL  HANK  ( RHO  / HI / M1F » H2 » H2P ) 

TYPE  243 / J f RHO  /HI/ HIT  <H2rH2P 

243  FORMAT ( 7 243  J / RHO / H 1 r HI F / H2 / H2P  7 /I OG  > 

Xi\TEH=RHO/X  ( N3  ) 

DO  123  I = 1 » NH2 

123  DUM  ( I » J > =AM  (I  t 21  H ) #XNS*H1+AM  <I/Jl+2>#<  XKTEM#XNS*HlF+Hl/<  2 . *XHS  > » 
TYPE  2 44 / ( ( BUM  £ I / J ) / J~ 1 / NM ) / I =1 / NM2 ) 

244  'FORMAT  ( 7 244  DUM  7 / ( 6G ) ) 


C 

C 

c 

c 


MULTIPLY  FROM  RIGHT  BY  COEFFICIENTS  FOR  INFINITE  SEGMENT  TO  DEFINE 
COEFFICIENTS  FOR  LAST  FINITE  SEGMENT 


127 


126 

543 

246 

247 


DO  126  I-l 7 NM 
I1--2*(I-1)  . 

3UM=0 . 

SUM2-0 . 

DO  127  J=1 f NM 

SUM  -SUM } BUM ( 1141 » J) *A ( J ) 

SUMS - S UH 2 4 DUM ( I H 2 / J > *A  < J ) 
AMI 1 1+1  * 1 ) -SUM 
AM \ 11+2  7 1 ) =SUM2 
ALFvMSr I )=SUM 
BET ( NS  r I ) ~SUM2 
TYPE  243  7-;  AM  ( 1. 1 ) /I  = l »NH2) 
FORMAT ( 7 245  AM7/120) 

TYPE  2467 (ALF (NS/ I ) » I=1/NM> 
FORMAT  < 7 246  ALFV6G) 

T YPE  247  » ( BET  ( NS  / 1 > / 1 -- 1 / NM ) 
FORMAT ( 7 247  EET  7 /6G ) 


L 

C 

L 

c 


FIND  COEFFICIENTS  OF  OTHER  SEGMENTS  BY  MATCHING  SCI 
USE  RELATION  IN  MIDDLE  OF  PAGE  43 


' 'ENTIALLY . . . 


NS2*NS-2 

IF (NS2 . LT. 1 > CO  TO  20 
IS=N3 

DO  123  IS2=1*NG2 

13=13-1 

ISM=IS-1 


u 

C 

C 

c 


MULTIFLY  FROM  LEFT  BY  SSRI'"' 
A I BOUNDARY  BETWEEN  CURRENT 


rcn  rr.Euicus  segment  evaluated 

:i  -PREVIC-US  OECMCNTS 


DO  12?  I -1 /NM 
II 72/11-1) 

Cl=C2C3(0Sf  IS»I)«<:E)  > 

S 1 • CS 1 : 1 < S T 1 1 3 » I > :'7?<  J I S ' > 

du  !;iik- ,1)  7d  :.*M7nn»: '+siham(ii+2»1> 

<29  DUM(lx+2.1)»QS<rs,  I)  t (-S'l*/4M(Ii4lii)+Ci+AM(IH2ll)  ) 


ft 

C. 

o 

o 

a 

e 

% 


0 

*> 

o 

o 

© 

© 


t I 


I 


*- 
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lOIS  PAGE  IS  BEST  QUALITY  PRACTICABLE 
COPY  PURIIISHED  lOM^Q  - 


TYPE  240>  (BUMCX  > 1 ) , I--1  >NM2) 
240  FORMA  M'  2-18  BUMV12G) 


MULTIPLY  FROM  LEFT  BY  U FOR  PREVIOUS  SEGMENT 

HU  130  I-lfNM 
U=2*<I-1> 

HO  130  J=1 »2 
SUM=0 . 

HO  131  K=1»NM 

131  CUM=SUM * US ( 1 3 > I » K ) *DUM  < 2* CK-1 ) f J > 1 ) 

130  AH ( 1 1 +2  > 1 ) =SUM 

TYPE  249, (AM( I > 1 > >1-1 , MM2) 

249  FORMAT  < ' 249  AM'/12G> 

MULTIPLY  FROM  LEFT  HY  U INVERSE  FOR  CURRENT  SEGMENT  • 

BO  132  1=1 >NH 
I1=2*(I-1> 

BO  132  J=l>2 
SUM=0 . 

BO  133  !\=1  >NM 

133  SUM=SUMH!IS<ISM, I ,K >*AM <2* <K-1 >+ J> 1 ) 

132  HUM< II >J» 1 )=SUM 

TYPE  230  > < BUM ( I > t ) > 1 = 1 > NM2 ) 

250  FORMAT ( ' 250  BUMV12G) 

OBTAIN  CGSFF.I  CENTS  FOR  CURRENT  SEGMENT  BY  MULTIPLYING  FROM 
LEFT  BY  SCRIFT-C  INVERSE  FOR  CUR'RENT'  SEGMENT  EVALUATED  AT 
RIGHT  OF  SEGMENT 


BO  134  I=-l  ,NM 
U=2*(I-1> 

ci=ccos ( os  < ism > j ) rx  < is> ) 

S1=C3IN < GS  < ISM , I > :::X ( IS)  ) 

SU>i=Cl*BUU < 11+1,1)  -SI #DUM (11+2,1)  /GS(  ISM>  I ) 
SUK2=31*DUM ( 1 1 + 1 > 1 ) +C1*BUM (1 1 +2  > 1 ) /OS ( ISM  > I ) 

AM (11+1,1; =SUM 
AM  < 1 1 +2  > 1 ) - SUM2 
. ALF(I3>r;=SUM  ■ 

134  BET •’  IS  > I ) =SUM2 

TYPE  251  > ( AM  ( I > 1 ) > 1= 1 > f.'M2 ) 

251  FORMAT  ('  251  AMV12G) 

TYPE  252 > < ALF < I S > I ) *I  = 1>NM) 

232  FORMAT ('  252  ALFV6G) 

' TYPE  253  » ( BET ( I S • I > » I =1 » NM ) 

253  FORMAT < ' 253  BETV6G) 

123  CONTINUE 

RECOMPUTE  COEFFICIENTS  IN  SOURCE  SEGMENT. 

TYPE  620. (QS(lfl) >1=1, NM) 

620  FORMAT ('  620  GS'/oG) 

BO  300  I =1 >NM 
11=2:)  ( l-l ) 

ci=ccos<Gsa>i>-:'X<i)) 

S1=CSIN<C3U.I>*X<1)> 

CUM  (1 1-M  1 1 )--Cl*AM  ( 1 1 + 1 > 1 ) tFt**M(  T l2,  t > 

600  BUK< I 1+2,1 >=03(1  >I)Y(-Sl-Vvat4’  . 1 ) »C1»AM<I1'?>1)) 

TYPE  601 » ; BUM  (I  > 1 ) » I = 1 . NM  •“  > 

601  FORMAT ('  601  DUMV125) 

BO  602  I ~ 1 » N M 

11  -2  MI-1) 

OC  602  J-l » 2 
SOMso. 

DO  603  ..... 


1 


ci  n n n n n 


3HIS  PAG®  IS  BEST  QUALITY  PRACTICABLE 

mil  ooire  pukiulshhd  ro  roc  — 


603  OOM -Clin »•-’«<  1 * X*K>*DUM<S1(<K-t  } hJt  1 > 

602  AlUIl  IJ,  J.)  -SUM 

C TYPE  601* (AM( 1*1)* I~1 *NM2> 

604  FORMAT < ' 604  AM'/12G) 

X25=SQRT(X( 1 ) ) 

LO  605  I--1  »NM 

RHO'-SCRT  ( XKM < I * 1 > >HX U > 

CALL  HANK <RH0*H1  *H1P*H2*!!2P) 

XK'TEii  -RHO/X  ( 1 ) 

Al=X2*rll 
A2^X2*H2 

A3--XKTCM  SX2*K1P  M !i/  ( 2 . *X2 ) 

A1 -XKTEM*X2*H2P »H2/(2. ?X2 ) 

A3-A1SA1-A2*A3 
TED -ril 
A1-A-4/A5 
A4-TEM/A5 
A2=-A2/A5 
A3=-A3/A5 

2UM (I1  + 1»1)*(AI #AM C 1 1 + 1 » 1 > +A2*AM < 1 1 +2 » ! ) ) / AO < I > 

605  HUM C 1 1 +2  * 1 ) - < A3  v AM i I H 1 * 1 ) +A4*AM  < 1 1 +2 . 1 ) > / AO  <1 > 

C TYPE  606*<DUiKI*l)*l-l*NM2> 

606  FORMAT < ' 606  DUM '/12G ) 

INPUT  RANGES  AND  DEPTHS  TOR  FINDING  PROP  LOSS 

20  TYPE  21 

21  FORMAT < ' INPUT  0 FOR  RANGE  VARIATION*  1 FOR  DEPTH ' VARIATION ' ) 
ACCEPT  2*IRZ 
IF ( IRZ.E3. 1 ) GO  TO  22 
TYPE  23 

23  FORMAT ! ' INITIAL  RANGE » FINAL  RANGE* ' * 

1'  STEP  SIZE') 

ACCEPT  2*R1*R2*BR 
Rl-Rl-itlOOO. 

RZ-RZ’MOOO . 

DR-DRKIOOO.  ■ 

R1--R1-DR 
F<-  R1 

MP?—  ( R2  -R ) /DR+ . 01 
NZ-1 

• DZ-O.  . 

GO  TO  24 

22  TYPE  25 

25  FORMAT < ' INPUT  RANGE*  INITIAL  DEPTH*  FINAL  DEPTH*  STEP  SIZE') 
ACCEPT  2*R*Z1»Z2*DZ 
Z1=Z1-DZ 

NZ=<Z2-Z1)/DZKC1 
NR=1 
DR -0 . 

24  DO  135  IR-1  * MR 
F;'R+DP 

ROUT -0/1009.  ' 

DC  136  IS-l.NS 

IFCR.LT.XCIG))  GO  TO  025  . ' 

136  CONTINUE 

RANGE  IS  IN  IfTIVTTE  SErpriir 


IS-NSP 

ZB1 - SCR  T <ZP<  t ) 1 

zi'N-snr.Tczr  cngd  > 

DO  137  I--1  *NM 

F.IIO  SORT < XKM  ( I * vsr ) ) *•' 

CALL  HANK{B.H0,Hl,HlPiM2iHlt) 


- i 1 - mis  FAQ £ IS  BEST  QUALITY  PRACTICABIl 

roou  copy  FUKfiisiiac  . ■ 


137  FR<I>--.MI)»II1/ZE: 

SUM=0. 

EC  750  1 = 1 -NM 

750  SUM -SUM ♦ KR  < I > *UR < I , NCP > 

SUiTCL'M/ZEN 

GO  TO  24 

S25  IF(IS.EO.l)  GO  TO  27 

RANGE  10  NOT  IN  EITHER  SOURCE  OR  INFINITE  SEGMENT 

ISM--IS-: 

X2=3uRT  <R)  . 

ZB 1 -SORT (ZE ( 1 ) > 

DO  13S  I :rl  ’NM 
GUM -0 • 

HO  13?  J~1  > NM 
GR=QS  < I CK » J > #R 

1 3?  GUM  =SL'f  W US  ( I SM , I , J ) * ( ALF  ( IS » J > *CC05  ( QR ) FEET  ( IS , J ) *C?  IN  (OR) ) /X2 
130  FR  ( I ) - EUK/ZD1 
GO  TO  26 

RANGE  IS  IN  SOURCE  SEGMENT 

27  ?B1=SQRT  (Z3 ( 1 ) ) 
no  140  1=1 »NM 
RHO-SGRT ( XKN (1,1)) KR 

CALL  HANK ( RHO , H 1 , !■'  1 F , H2 , 1 12P  > 

140  FR  C I ) =A0 (I ) * ( ALF (1,1) Kill  FEET (1,1) *K2 ) /ZE1 

OUTPUT  NODE  INFORMATION 

• • 

26  IF(NOUT.Ea.O)  GO  TO  45 
J=1 

HG  150  1=1, NM 

IF ( I . ME . IOUT C J ) 5 GO  TO  }50 

I1=2*(J-1) 

G< 11+1 )=CABS(FR( I ) ) 

G( I1F2 ) =CLOG (FR ( I ) /G ( 11+1 ) ) /EYE* 130 . /PI 
IF(J.EG.NOUT)  GO  TO  46 
J= JF1 

150  CONTINUE 

46  WRITE  (23,47)  ROUT , ( G (I ) , 1=1 , N0UT2 ) 

47  FORMAT < FI 0 . 3 , 1P6E 10 . 3/ < 10X , 1P6E1 0 . 3 ) ) 

CONFUTE  PROP  LOSS 

45  IF  (IS.  EC!  .NS?)  CO  TO  751 
ISP=I3!1 

FRAC= ( P-XNAG ( IS ) ) / ( XNAG (ISP) -XN  AG (IG ' ' 

SUM=0 . 

ZK1=SGRT(IB(IS) ) 

ZB2=SC1RT (ZB(  ISP)  > 
no  142  1=1, NM 

142  SUM  ■ SUM  : FR  ( I ) * ( UR  ( I , I S ) /ZBH  FRACf  ( UR  N-FU.  '702-UR  ( T , IS  ) /ZB1 ) ) 

751  TL=-20.*AL0C10 ( 4 . KPIKCAES ( SUM ) ) 

IF  ( IRZ.EC!  • 1 ) GO  TO  2C 

IF  (ROUT  . NE . 1 . .OP . ROUT . ME.  10  . . OF  • • • 75 . . OR. ROUT . ME. 30.  .OR  . T 

1.NE.75.  .Oii. ROUT. NE.  100. ) GO  T3  ' 

TTFE  29, ROUT, TL 
2?  FORMAT (TIC. 5, FI 0.2) 

1234  CONTINUE 

WRITE  (22,???)  PCUT , TL 
V 7?  FORMAT <ri3.S,FlC.«> 

Cij  1 U 123 

28  rvrr  -*?  - .:r  ,t». 

WRITE  ( Z2, 29  ) ZD,TL 


■"3. .OR. ROUT. ME. 30. .OR. ROUT 
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IBIS  PAGE  IS  BEST  QUALITY  PRACTICABLE 
PROtf  COPY  FURBISHED  TO  LD£ 


us>  Cunhhuu  — * — ~~t-  -•  •-— — • - ~ 

T“-l . 

write  t22.29>  t 
TYPE  30 

30  F03tt.YT('  INPUT  0 FCR  NEW  RANGE  GR  PERTH > 1 FOR  END') 
ACCEPT  2.IEND 
IF(IENO.EC.O)  GO  TO  20 
STOP 
END 


1 1 


THIS  PAGE  IS  BEST  QUALITY  PRACTICABLE 
FROM  COPY  FURNISHED  TO  DOC 


SUBROUTINE  NAGL < NM » NLP ,NSP» NMODE >2»ZB.Xi;, XKNS t UOS . UPS , 2S » XKi ! • UO . 
luh  ) 

DIMENSION  2 < •I ) .ZB<4)  » XK<4.4>  »XKN(  15.4)  *IJO(  13)  ,UR(15.4>  »NM0DC(4) 
DIMENSION  XNN3  ( IS  » 4 ) . UOS  < 15 ) »Uf\S  ( 15. 4 ) 

DOUBLE  PRECISION  UO.r  UP  > XKN 
DO  100  1SP-1.NSP 
IM=KUODE(ISP> 

IF(IH.ECl.O)  GO  TO  GOO 

IF < ISP. ECJ . 1 )RE AD< 23 . 1 > <U0< I > » 1=1 » IM) 

1 FORMAT < 10CD) 

C IP ( ISP. EG  . 1 )TYPE  237 , (UO(I).I-l.IM) 

237  FORMAT ( ' 237  U0'/<3G>> 

READ <27. 1 > ( UR < I . ISP > » 1 = 1 , IM > 

C TYPE  238.  (UR( I » ISP) »I=1 > IM) 

238  FORMAT ('  23S  URV3G) 

READ ( 23 , 1 X XKN ( I , ISP ) » 1 = 1 , IM  > 

DO  101  1=1. IM 

101  XKN ( I » ISP) =XKN < I , I SP ) **2 
C TYPE  239.<XKN(I»ISP).I=1,IH> 

239  FORMAT ('  239  XKNV3G) 

IFtZS.GT.O. .OR.ISP.NE.l)  GO  TO  100 
DO  302  II.  IM 

302  UO<I)=UR<Ifl> 

GO  TO  100 

300  IM'-NMODE  ( i ) 

NMODE ( ISP) “NMODE ( 1 ) 

DO  3C1  1=1. IM 
UR ( 1 . 1 5P  > =UR ( 1 . 1 > 

301  XKN ( I . I 3P ) =XKN ( I . 1 ) 

100  CONTINUE 

DO  110  ISP= 1 » NSP 
IM=NMODE ( ISP ) 

DO  110  1=1. IM 
UOS(I)=UD(I) 

URS< I » ISP  >=UR( I » ISP) 

110  XKNS( I » ISP)=XKN ( I » ISP)  ' 

RETURN 


END 


nn on  nnnnn  nnnnn  non  nno 


IBIS  PAGE  IS  BEST  QUALITY  FRACIICABI* 

• - '•  FROM  COFY  FURNISHED  10  DOC — 

, TY  FR\R\"C 

* I 

.TTY  FORM 
I .TY  FORM 

SUBROUTINE  FORM ( IS . U. UI . 0 . ICOUP.NMODE > 

' COMMON  /DEPTH/  XNN<15.4)» 

l iuo(is> 

« COMMON  /FORMC/  X ( 3 ) . Z < 4 ) »ZD ( 4 ) » NM . ML . NLM .NLP . XK ( 4 . 4 ) . DEN< 4 ) 

1 t XNAG ( 4 ) 

DIMENSION  AC (3) »ST<13> .CT(13> .C(15> 10) , V( 15. 15) . DUM< 15. 15> . 
lS<15.15).SI(15»15)»EIG(iO> .U<15, 15  > , UI  < IS / IS)  .0(13)  »DHAT<  15. 13)  t 
2EHAT  ( 120). DVEC(  15. 13).  DEIG<  15) , UK  (-13).  UK2(1480)  » 510(15,15) . 

3XMXEL (15.15) . NrtODZ  < 4 ) 

COMPLEX  SUM . C . V . BUM > 5 . SI » EIG . U . UI . « . DM AT . El lAT . DVEC . 3 1 G 
DOUBLE  PRECISION  XNAG 

COMPUTES  U=S*SIGMA*S1  AND  ITS  INVERSE  AND  THE  EIGENVALUES  Q OF 
FINAL  TRAN3F  QZMAT-ION  . 

MATSIZ=15 
NP=KM*(NM+l)/2 
NQ=3HNM 

NUK2=2*NM*(NMH> 

ISM=IS~1 
ISP=IS+1 

LINEARIZE  KN*#2  IN  RANGE 
DO  100  1=1. NM 

ST ( I )=(XKN ( I . ISP)-XKN( I.IS))/(X(IS)-X(ISM) ) 

100  CT ( I )=XKN( I » IS) -ST (I)*X(I5M) 

TYPE  202. C3T(I) .1 -1 .NM) 

202  FORMAT ( ' 202-  ST'/3C> 

TYPE  203. <CT( I ) . 1=1 .NM) 

203  FORMAT ( ' 203  CT'/3G) 

LINEARIZE  K**2  IN  DEPTH  AND  RANGE  • 

DO  101  IL=1.NLM 

101  AC(IL)=(XK(IL+1» ISP  >-XK  (•IL.ISP)-(XK(IL+1,I3)-XK(IL»IS)))/ 

1(  <X(IS)-»X(ISM)  >*(Z(IL+1)-ZCIL) > ) 

AC(NL)  = ( <XI\  (NLP. ISF ) -XK(NL. ISP ) )/(ZB( ISP)-Z(NL) ) 

1-<XN(NLP»  IS  )-XK  (NL » IS ) >/.<ZB  ( 15 ) -Z(NL  ) ) )/<X(I3)-X(ISM)  ) 

. TYPE  204. ( AC ( I ) . 1=1 »NL ) 

204  FORMAT ('  204  ACV3G) 

FORM  MATRIX  C ON  PAGE  37 

CALL  MXEL ( NMODE , IS . XNAG . XMXCL . ZD ) 

TYPE  203. < (XMXEL ( I . J) . J"1 .NM) .1=1 »MM) 

205  FORMAT ( ' 203  XMXEL' /< 30)) 

DO  104  1=2. NM 
IM=I-1 

DO  102  J=1»IM 

C ( I . J )=-XMXEL ( I . J ) / ( (ST ( I ) -ST ( J) >*<X<I3)+X(ISM> )/2.+CT(I)-CT< J) ) 

102  C<J.I)=-C(I.J) 

104  C(I.I)-0. 

C(1»1)=0. 

TYF  E 204.  < < C ( I * J ) . J - 1 .NM  > . I 1 / NM  ) 

204  FORMAT ( ' 204  C'/(40>> 

FORM  MATRIX  V ON  PACE  34 


DO  105  1=1 .NM 
DO  103  J=1»I 

SL’M=0. 


-?>5*- 


__  *|£T  QjjAim  moneABii 

SSSJ-SU®®'1 — ' 


HO  107  K-l , NM 

IFCK.EO.I.OR.K.EO. J>  GO  TO  107 
SUM-SUM  t-C  < I » K ) *C  ( K , J > 

107  CONTI NUZ 
SUM=-SUM 

IF(I.ECJ.J)  GO  TO  2 

V(I» J)=-C(I, J)*(ST(I)  -ST(J))/((ST(I)-ST<J)  >*(X(IS)+X(ISM) )/2. 
1+CT(I)-CT< J) > +SUM 

U(J,I)=-C(  J,I)*(ST<  J)-ST(I)  )/(  <ST(  J)-ST(I)  )*<X(IS>  f-X(ISM)  )/2. 
1+CT  < J > -CT  < I ) ) +SUM 
GO  TO  105 
2 V ( I , I ) —SUM 
105  CONTINUE 

IF ( ICOUF .EO. 1 ) GO  TO  302 
HO  803  1=1 »NM 
no  804  J=1 » NM 
U(I, J>=0. 

U(Jr I)=0. 

UI(I,J>=0. 

804  UI(J,I)=0. 

U'.M)»1. 

UI(Z>I)-1. 

a(I)=SGRT<  < XKN (I»IS) +XKM < I * ISP ) ) /2 » ) 

803  CONTINUE 

TYPE  803, ( (U(I, J) *J=1»NM) ,I=1,NM) 

805  FORMAT ('  305  U'/<60>> 

TYFE  806,<(UI(I,J>,J=1»NM),I*1,NM> 

806  FORMAT  ( ' 806  UIVC6G)) 

TYFE  307, (0(1) ,1=1, NM) 

807  FORMAT ('  007  Q'/6p) 

RETURN 

802  CONTINUE 

TYPE  207 ,((V(I,J), J=1 , NM ) , 1=1 , NM) 

207  FORMAT <'  207  V'/C6G)> 

V 

DIAGONALIZE  C. 

C IS  ANTISYMMETRIC  (ANTIHERMITIAN) • . .EIGENVALUES  ARE  IMAGINARY 

CALL  HMD  I AG  ( NM , NP , NO,  C , EIG , 3 , SI , OMAT , EMAT , DVCC , DEIG,  UK ) 

TYPE  208, ( (SCI, J) ,J=1,NM),I=1,NM) 

208  FORMAT ('  203  S'/(6G)) 

TYPE  209,  ( (SKI,  J)  ,J=1»NM)  ,I  = 1,NM) 

209  FORMAT ('  209  SI'/(6G)) 

• TYPE  210, (EIG(I), 1=1, NM) 

210  FORMAT ('  210  EIGV6G) 

DEFINE  DIAGONAL  MATRIX  SIGMA 
DO  108  1=1 , NM 

SIG ( I r I > --CEXP (-EIG(I)T(X(IS)+X(ISM) )/2. ) 

IrU.EO.l)  CO  TO  103 
IM--I-1 

DO  300  J-’  1 * I M 
SIG ( I , J) =0 . 

300  SIG ( J , I ) - 0 . 

108  CONTINUE 

TYPE  211 , ( (SIGCI ,J> ,J'1,NM) ,I=1,NM> 

211  FORMAT ( ' 211  SIGV(6G>) 

FORM  MATRIX  UfK 
DO  10?  1*1 » NM 

109  U(I,I)-VU ,1) ! ( XKN ( I , IS ) fXKN(I.ICr>>/2. 

. Tttr  21Zr(<V( i, j),j  i»mn),t  i,nh) 


112  FORMAT  (*  H2  \J'/(66)) 


/ 
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DO  HO  1 = 1 fNM 
DO  110  J=1 pNM 
GUM-0 . 1 
DO  111  K=1fNM 

111  SUM -SUM  fV  (I » K ) :*S  ( K f J ) 

110  DUIK  I f J ) =SUM 

TYPE  213f  < <DUM(IfJ),J=1fNM>fI=1fNM) 

213  FORMAT  ('  213  DUM'/(6G)) 

FORM  MATRIX  T ( CALL  IT  0 TO  SAVE  SPACE)  ON  PAGE  38 

DO  114  1=1 fNH 
DO  112  J=1/NU 
SUM=0. 

DO  113  K=1fNM 

113  SUM  =SUMFSI < I fK ) TDUMCK f J > 

112  V(If J)=SUM*CEXP< (EIG(I)-EIG(J) >#<X( IS)+X( ISM) )/2. ) 

114  V<IfI>=V<I»1)-EICCI>#*2 

TYPE  21 4 f < ( V(I f J) f J=1 fNM) f 1=1 fNM) 

214  FORMAT ( ' 214  V'/<83)) 

DIAGONALIZE  T ( CALLED  V IN  CODE) 

CALL  EIGCC(VfNMfMATSIZf2fCJfUfMATSIZfWK2fIER> 

TYPE  £88 f (0 < I ) f 1=1 fNM) 

888  FORMAT ( ' 833  C!'/6G) 

DO  450  1=1 fNM 
450  Q ( I ) '-CSQRT  (0(D) 

TYPE  215f (0(1) fI=1fNM) 

215  FORMAT  ('  215  C1'/6G) 

TYPE  216 f ( ( U ( I f J ) f J=1fNM>  fI  = 1fNM) 

216  FORMAT < ' 216  U'/(6B)> 

THE  FOLLOWING  CARD  IS  TEMPORARY  TO  CHECK  PERFORMANCE  OF  EIGCC 
TYPE  400  f WK2 ( 1 ) 

IF(UK2<1) .GT.i . > TYPE  400fUK2(1) 

400  FORMAT < ' EIGCC  PERFORMANCES r 1PE10.3) 

FIND  INVERSE  OF  U 

DO  410  1=1 fNM 
DO  410  J=1 f NM 
410  V ( I f J ) =U ( I f J ) 

CALL  LEQTOCCVfNMfKATSIZfDUMfIfMATSIZfIfUKZfIER) 

IFCIER.NE.129)  GO  TO  412 
TYPE  411 

411 • FORMAT < ' U IS  SINGULAR') 

STOP 

412  IF ( IER.EO . 130 > TYPE  413 

413  FORMAT ('  ITERATION  FAILED  TO  IMPROVE  U INVERSE.  MATRIX  IS'f 
1'  TOO  ILL  CONDITIONED') 

DO  414  IM=1fNM 
DO  415  JM=i fNM 

415  DUM  < JM  f 1 ) =0 . 

DUM(IMf1)  = 1 . 

CALL  LE0T2C < V » NM . MATS IZ f DUM f 1 f MATS I Zf2f WK2 f I ER ) 

IF< IER.NE. 12/)  GO  TO  416 

TYPE  411 

STOP 

416  IF < IER.EO. 130)  TYPE  413  * 

DO  417  JM  -1 f NM 

417  UKJMfIM)  DUM(JMfI) 

414  COPT  r NIT-  _ 
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C TYPE  217, <CUI<I,J.',J=1,NM>,I  = 1»NM> 

C ' 217  FORMAT ( * 217  UI'/(6G)> 

FORM  MATRIX  SIGMA*S1 
HO  115  1=1 »NM 
no  115  J=1 »NM 

115  nUM<I»J)=3IG(If I>*UCI»J> 

TYPE  218,  < (DUM<I , J) , J=1 ,NM),I=1,NM> 

218  FORMAT C ' 218  DUM'/<6G)> 

FORM  MATRIX  U=S*SIGMA*S1 

DO  117* 1=1 » MM 
HO  117  J=1 , NM 
SUM =0 . 

no  118  K-l ,NM 

118  SUM  -SUMPS < I »N> *iiUM (N»  J> 

117  U ( I , J ) -SUM 

TYPE  219, ( <U( I » J) , J=1,MM> ,1=1 ,NM) 

219  FORMAT ('  219  U'/(6G>> 

FORM  MATRIX  SI  INVERSE*SIGHA  INVERSE 

•DO  119  1=1, NM 
DO  11?  J=1 , NM 

119  DUMd,  J)=UI(I,J)/SIG<J,J) 

TYPE  220, ( (DUM (I»J),J=1»NM)»I=1»NM) 

220  FORMAT  ('  220  DUM'/dO" 

FORM  MATRIX  UI  INVERSE=S1  INVERSE*SIGMA  INVERSE*S  INVERSE 

• • 

DO  120  1 = 1,  NM 

i DO  120  J=1 » NM 

\ | SUM=0. 

j DO  121  K=1 , NM 

I 121  SUM=SUM+nU:i<I,K>*SI<K,J> 

120  UI<I,J)=S’JM 

T C TYPE  221 , < (UI (I, J) > J-l ,NM) , 1=1 »NM> 

221  FORMAT ('  221  UI'/<6G>> 

RETURN 

. END 


■ c 


IHIS  PAGE  IS  BEST  QUALITY  PRACTICABLE 
FROU  OOP!  FURNISHED  TO  DDC ' 

TY  COEF1 

SUBRCUT I NE  COEF 1 ( NM  , AM  r V , A , AO , ALF » PET  » FRAC) 

DIMENSION  AM (30. 30) .0(15) , A (IS) .ACUCI) »ALF<3, 15) , DET (3,15) , 
1C<15,15)»V1(15),V2(15)»V3(15) 

COMPLEX  AM » V » A » AO , ALF » BET » C , SUM  , SUM2 

ROUTINE  COMPUTES  COEFFICIENTS  IN  INFINITE  AND  SOURCE  SEGMENTS  BY 
MATCHING*  USING  AN  ITERATIVE  APPROACH. 


NC=10 
TYPE  911 

911  FORMAT ( ' INPUT  MAXIMUM  NUMBER  OF  ITERATIONS') 

ACCEPT  912. NC 

912  FORMAT (G) 

FOR  MODES  WHICH  DO  NOT  PROPAGATE  AT  SOURCE.  V=0. 

SET  V- • 00001 ^SMALLEST  NONZERO  V FOR  ITERATION  PROCESS. 

VMIN=CABS(V(1)) 

DO  100  1-2, NM 
VT=CA3S(V(I)) 

IF(VT.LT.VMIN.AND. VT.GT.O. ) VMIN=VT 

100  CONTINUE 
VMIN=VNIN#l.E-5 
DO  101  I-l.NM 

IF(V(I).EQ.(0.»0.))  V(I)=VMIN/(0. ,1. ) 

V 1 ( I ) =RE AL  < V ( I > > 

V2( I )=AIMAG( V( I > ) 

Vo  ( I > "FRACSCABS < V (I ) ) 

101  CONTINUE 

TYPE  229. (V(I>, I-l.NM) 

229  FORMAT  ('  22?  V'/3G) 

TYPE  230. (V1(J)»I=1 »NM) 

230  FORMAT ( ' 230  VI' /3G) 

TYPE  231 . (V2(I) .1=1 >NM) 

231  FORMAT < ' 231  V2'/3G) 

TYPE  232. (V3( I > ,1=1 »NM) 

232- FORMAT  ('  232  V3'/3G> 

COMBINE  ELEMENTS  OF  AM  FOR  USE  IN  EQUATION  4 OF  PAGE  NAGL.4 

DO  102  1=1, NM 
I1=2*(I-1) 

DO  102  J=1 , NM 

102  C(I»J)=AM(I1+1,J)-AM(I1+2,J> 

TYPE  233, ( <C(I, J) , J=1,NM) ,I=1»NM) 

233  FORMAT ('  233  C'/(4G>) 

DEFINE  A FOR  ZEROTH  ORDER  (PAGE  NAGL.5) 


£>  i 


DO  103  1=1, NM 
103  A( I ) =V ( I ) /C  < I , I ) 

TYPE  234, (A (I), 1=1. NM ) 

234  FORMAT ( ' 234  A'/00) 

IC-ITERATION  COUNTER 

IC=1 

CO  TO  SO 

DEFINE  A FOR  IC-TH  ITERATION 

4 DO  10*.  1 = 1 »NM 
SUM- 0. 

DO  100  J 1 »NM 

Q,S)  GO  TO  l'"'S 
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SUM-SUM  ) C ( I , J ) #A ( J ) 

105  CONTINUE 

104  A(I)=(V(I) -SUM ) /C ( I , I ) 

THE  FOLLOWING  CARD  IS  TEMPORARY  TO  TEST  ITERATION  PROCEDURE 

50  IW=0 

DO  106  1=1 »NM 
SUM=0. 

DO  107  J=1 , NM 
107  SUM=SUM+C<I, J)*A(J> 

X1=REAL < SUM) 

X2=AIMAG  < SUM ) 

THE  FOLLOWING  2 CARDS  SHOULD  BE  INCLUDED  WHEN  ITEERATION  TEST 
13  COMPLETED 

IF ( ABS ( Xl-Vl (I)).GT.V3(I))  GO  TO  1 
IF< ABS ( X2-V2 (I)).GT.V3(I>>  GO  TO  1 

THE  FOLLOWING  3 CARDS  ARE  TEMPORARY  TO  TEST  ITERATION  PROCEDURE 

IF( ABS < Xl-Vl < I ) ) i 3T . V3< I > ) IU=1 
IF  < ABS<X2-V2(  I)).GT.V3(I))  IW=1 
TYPE  875 , lCrlfA(I) , XI » X2, VI < I ) , V2 ( I ) 

875  FORMAT  <213 , 1P2E15 . 7/6X , 1P4E15 .7  > 

106  CONTINUE 

THE  FOLLOWING  CARD  SHOULD  BE  INCLUDED  WHEN  ITERATION  TEST  IS 
COMPLETED 

GO  TO  2 

THE  FOLLOWING  CARD  IS  TEMPORARY  TO  TEST  ITERATION  PROCEDURE 

IF<IW.EQ.O>  GO  TO  2 , 

1 IF < IC.EQ.NC)  GO  TO  3 

START  NEW  ITERATION 

• IC=IC+1 
GO  TO  4 

MAXIMUM  NUMBER  OF  ITERATIONS  REACHED  WITHOUT  CONVERGING 
3 TYPE  5 

5 FORMAT (//'  A NOT  CONVERGED.  SUM,  V='> 

DO  109  1=1 ,NM 
SUM=0. 

DO  110  J=i ,NM 
110  SUM=SUM+C<If J)*A< J) 

TYPE  375 ,IfICfA(I)f SUM , V(  I ) 

109  CONTINUE 
TYPE  600 

600  FORMAT ( ' TYPE . . . -1  STOP'/ 

1'  0 USE  LAST  VALUES  FOR  COEFFICIENTS'/ 

2'  N ITERATE  N MORE  TIMES') 

ACCEPT  601 fNC 

601  FORMAT (U) 

IF (NC)  -602,2,603 

602  STOP 

603  IC-1 

GO  TO  4 

; ist Mr  r.ouvrpocn  vm  ur.s  of  a <or  vahts  from  mc-th  iteration), 

DEFINE  AOrAV<i  r.l.l  , M.  M FCR  • ’•  ' .i'  I fiT.V.M 
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TY 

SUBROUT  I NE  HMD I AG  < N * Nr » NCI  t CHAT  r CF I G » CVCC » C VECC » DM  AT  » EMAT » DVEC , 

1 DE IG > UK ) 

COMPLEX  CM AT  » CEIG  * CVEC » CVECC  , E YE 
COMPLEX  DMAT. EMAT. DVEC 

DIMENSION  CMAT  ( 15.  IE)  ?CVEC<  15*  13)  » CVECC < 15, 15)  .CEIGUS) 

DIMENSION  DMAT (N.N ) .EMAT < NP ) » DVEC ( N.N) »DEIG(N) , UK(NO) 

MUST  DIVIDE  BY  EYE  BECAUSE  MATRIX  IS  ANTISYMMETRIC. 

THIS  MAKES  MATRIX  HERMITIAN  UHICH  CAN  BE  EASILY 

DIAGONALIZED.  ORIGINAL  MATRIX  HAS  EIGENVALUES  WITH  SAME  MAGNITUDE 
AS  THE  HERMITIAN  MATRIX  BUT  ITS  EIGENVALUES  ARE  PURELY  IMAGINARY. 

THE  EIGENVECTORS  ARE  THE  SAME. 

EYE=(0 • * 1 . ) 

DO  150  I = 1 » N 
DO  150  J~1 j N 

150  DMAT ( I , J)  -GMAT ( I , J) /EYE 
CALL  VCVTCH ( DMAT , N.N, EMAT ) 

CALL  EIGCH (EMAT »N»  1 .DEIG.DVEC . N, WK» IER) 

MUST  MULTIPLY  BY  EYE  BECAUSE  ANTISYMMETRIC  MATRIX  HAS  IMAGINARY 
EIGENVALUES 

DO  151  I-l.N 
CEIG(I)=DEIG(I)*EYE 
DO  151  J=1 » M 
CVEC  ( I , J)=IiVEC  < I » J) 

151  C VECC  < J » I > =CON JG  < CVEC ( I , J > ) 

RETURN 
END 

«■*  i - 
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TY  MXCL 
00050 
00100 
00200 
15) 

00250 

h0DE<4) 
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SUPROUT  TNE  MXEL  ( .VNODE  * NPRS  * D * OXEL , OZB ) 

IMPLICIT  DOUBLE  PRECISION  (A-ll-p-Z) 

DIMENSION  X( 15* 1) *A( 15*4) *B<15*4>  *U(15*4>  *UP(15* 4) *UINT( 
DIMENSION  D(4) *SL<4>  »AL(4)*S(15*15> ,0XEL(15*15) *0ZD<4>*N 


t 

00250 

M-HMODE  ( NPRS  I NC-- 1 ) 

%• 

00350 

DO  36  MC=1 *2 

00360 

DO  35  NP=1*NPRS+NC-1 

C 

00300 

DO  30  J-l,M 

00400 

READ (21 *40)  NB 

00500 

DO  30  1-1 *NB 

C 

00600 

30 

READ <21*10)  X( J*I).A( J*I>  *B< J,I> ,U< J*I>  *UP( J*I) 

00700 

10 

FORMAT (SD) 

00800 

40 

FORMAT (14) 

c 

00900 

C 

TYPE  80 *( (U(J*I)*I=1*4)*J=1*M) 

01000 

80 

FORMAT (4D16. 7) 

01100 

C 

TYPE  80* ( <UP<J,I) *1=1*4) *J=1*M> 

• t 

01200 

35 

CONTINUE 

01300 

IF(NC.ECJ.2)  GO  TO  39 

01400 

READ (21 *40)  NB1 

A 

c 

01500 

DO  31  I==1»NB1 

01600 

31 

READ  < 3*1  *10)  Cl  * AL  ( I ) * C2  *C3  *C4 

4 

01700 

DO  33  i\L=2*NB 

1 

c 

' 01800 

38 

SL ( KL ) = ( AL ( KL ) **3-A ( 1 * KL ) **3 ) / 

01900 

1 

< D < NPRS  H ) -D ( NPRS ) ) 

* 

02000 

39 

CONTINUE 

c 

02100 

CLOSE (UNIT=21) 

f 

02600 

DO  100  J=2*M 

02700 

DO  100  K=1  * J--1 

c 

02800 

IFCNC.EO.l)  S ( J * K > =0 • 

02900 

DO  110  1=2  * NB 

(i 

03000 

BP=(B< J*I)-B<K*I) >*A(JfI>**2 

c 

03100 

BS-B<  J* I )4B<K* I ) 

* 

03200 

V=A  ( J * I ) * ( BS+2  . *X  ( J,  I ) ) *U  ( J,  I > *U  ( K » I ) +2  . *UP  ( J*  I ) *IJP  < K * I > 

\ 

03300 

U-(BD#X< J*I)-2.#A( J*I>/BD>*(UP(J*I)*U(K*I)-U(J*I>*UP(Kf I 

c 

) ) 

03400. 

Y=A<J*I>*(BS+2.*X(.J*I-1)>*U(J*I-1)#U(K*I-1)+ 

03500 

1 

2.*UP<  J*I-1)*UP(K*I-1>#<A<  J*I-1)/A(  J»I)  )**2 

i 

4 

c 

03600 

Z=(X(J*I-1  )#BD-2.  :!<A(  J*  I > /BO  >#  (UP  < J*  1-1 ) *U<  K*  1-1 )- 

1 

03700 

1 

U(  J*I-1  )*UF(K,  1-1 ) )*A<  J*  1-1 ) / A ( J * I ) 

I ' 

03800 

S('J»K)=S(  J*K)  + ( V+W-Y-Z)/BB/JSD')<SL<  I )/2./(0ZB(NPRS+NC-l ) )# 

9 

c 

*2* 

% 

, 

03850 

OXEL( J»K)=S< J*K> 

03860 

105 

FORMAT (214 *D) 

, 

c 

03870 

110 

CONTINUE 

* 

039 00 

100 

CONTINUE 

4 

03950 

36 

CONTINUE 

c 

04000 

CONTINUE 

• 

02500 

RETURN 

* 

02 600  ’ 

END 

* 

c 

• 

' 

c 

c 

c 

. 

I 1 

o 
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TY  HANK 

SUBROUTINE  HANK < RHO . HI » H 1 P , H2'» H2P > 

COMPLEX  HI »H1P»H2»H2P 

CALL  HANLO < RHO » HI » HIP » H2 » H2P ) 

RETURN 

END 


* 


j t HANLO 
OOIOO 
00200 
00300 
00400 
00500 
00600 
00700 
00800 
00900 
01000 
01100 
01200 
01300 
01400 
01500 
01600  SC 
01700 
01800 
01900 
02000 
02100 
02200 
02300 


t H2P ) 


SUBROUTINE  HANLG C X » HI . HIP » H2 
REAL  JOfJl 

DOUBLE  PRECISION  PI»2»A»P»Q»R»SrX 
COMPLEX  HI fH1P»H2»K2P 
PI=3. 1415726535? 

Z=3*X 

Z1=DSCJRT  < 2/PI/X ) 

A=X-PI/4 

P=l-4 . 5/Z+ #2+9*25*49/24 . /Z**4 

Q=-l/Z+9*25/6 . /Z**3 

R=l+7  « 5/Z#*2 

S=3/Z-9*35/6 . /Z **3 

J0=Z1*SNGL  < F*DC05 ( A) ~0*DSIN  *A ) ) 

YO=Z  1 *SNGL  ( P*BSI N ( A ) +Q*DCOS  ( A ) ) 

Hl=J0+<0. >1. >*Y0 

FORMAT <2G> 

J0=-Z1*  < R#DSIN<  A ) + 3 *DCOS( A) ) 
YO=Zl#<r<#DCOSC  A)-S’?DSIN(  A) > 
HlP=J0+<0. »1. >*Y0 
H2=C0NJG(H1) 

H2P=C0NJG<H1P> 

RETURN 

END 


5 

3> 
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TY  DAIRY 

p 

00010 

00020 

9 

00030 

C 

IRY 

00040 

c 

9 

00050 

c 

9 

00060 

c 

00070 

c 

9 

00080 

c 

9 

00090 

c 

00100 

c 

9 

00110 

9 

00120 

c 

00130 

c 

# 

00140 

c 

e 

00150 

c 

00160 

c 

*s 

00170 

c 

f 

00180 

c 

00190 

c 

9 

00200 

c 

c 

00210 

c 

00220 

c 

0 

00230 

c 

0 

00240 

00250 

0 

00260 

0 

00270 

00280 

0 

00290 

0 

00300 

c 

00310 

0 

00320 

4L 

00330 

• 

* 
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SUBROUTINE  DAIRY < DX f AI fAIPfBI f DIP) 

26230 

IMPLICIT  DOUBLE  PRECISION  <A-HfO-Z> 

FOR  DOUBLE  PRECISION  ARGUMENTS r THIS  ROUTINE  CALCULATES  THE  A 
26240 

FUNCTION  AI (X)  AND  ITS  DERIVATIVE  AIP(X).  IT  ALSO  FINDS 
26250 

THE  OTHER  REAL  LINEARLY  INDEPENDENT  SOLUTION  BI<X)  AND 
26260 

ITS  DERIVATIVE  BIP<X>. 

26270 

THE  DEFINITIONS  AND  NORMALIZATIONS  ARE  AS  IN  NDS  HANDBOOK 
26230 

OF  MATHEMATICAL  FUNCTIONS f P . 446 
26290 

THE  METHODS  USED  ARE  POWER  SERIES  EXPANSION  FOR  SMALL  X 
26300 

AND  GAUSSIAN  INTEGRATION  FOR  LARGE  X 
26310 

DIMENSION  X<16)fW<16)fXSQ(16) 

26320 

DOUBLE  PRECISION  DXf AI fAIPfBI fBIP 
26330 

DOUBLE  PRECISION  XS  fXCUBEfAISUMf AIPSUM 
26340 

DOUBLE  PRECISION  DFfDFPfDGfDGP 
26350 

DOUBLE  PRECISION  FJM2f FJM1 fFJ fFJPI f FJP2f FACTOR 
26360 

DOUBLE  PRECISION  C1fC2fR00T3 
26370 

DOUBLE  PRECISION  DZETAfDARGfDROOTX 
26380 

DOUBLE  PRECISION  R00T4X f S f CO f RATIO f EFAC f ZETASQ 
26390  v 

DOUBLE  PRECISION  SUMRfSUMIfSUMRPfSUMIPfTERMRfTERMI 
26400 

DOUBLE  PRECISION  DZEROfDAfDBfDENfONE 
26410 

DOUBLE  PRECISION  XfUfXSQ 
26420 

DOUBLE  PRECISION  RSOfTEMPf  RTPIfRTPI2 
26430 

DOUBLE  PRECISION  TERMAf TERMS 
26440 

LOGICAL  NEEDS I 
26450 

DATA  DZEROfONE  /O . ODOf 1 .ODO/ 

26460 

DATA  R00T3/1 .73SCGnao7G63077DO/ 

26470 

DATA  Cl fC2  /.  3550230538370  1/T'Of  .253S19403792807D0/ 

26480 

DATA  RTPI  /. 20209  *.771773373  IDO/  ' 

26490 

DATA  RTPI2/.564189533-..}77S62r'0/ 

26500 

POSITIONS  AND  WEIGHTS  FOR  10-TERM  SUM  FOR  AIRY  FUNCTIONS 
26510 

DATA  U<  1)  / 3.  lS-lS"’ : ■' 7.1276 4 787D-  14/ 

26520 

DATA  W<  2)  / 6.639-*::-  '1753492111-11/ 

26530 

DATA  U(  3)  / 1 . 7503:*  * . t •* ’.5669D  03/ 

26540 
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1 

00480 

00490 

C, 

\ 

00500  ' 

'»  ; 

00510 

00520 

<0 

! 

00530 

<0  j 

00540 

00550 

(0 

00560 

( 0 ! 

00570 

j • 

00580 

<t 

00590 

(0 

00600 

00610 

(0  ‘ 

00620 

< 0 

00630 

00640 

(0 

00650 

' * 
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DATA  U<  4)  / 1.3VirS?23?04353l8D~06/ 

26550 

DATA  UK  5)  / 4.4.1S0966639204350D -05/ 

26560 

DATA  UK  6)  / 7* 1535010917718255D-04/ 

26570 

DATA  W<  7)  / 6.49593661033353010-03/ 

26580 

DATA  UK  9)  / 3.64404153757732220-02/ 

26590 

DATA  UK  9)  / 1.43997924105909990-01/ 

26600 

DATA  U<10>  / 8.12311413362614360-01/ 

26610- 

DATA  X<  1)  / 1.40330810721809640+01/ 

26620 

DATA  X<  2)  / 1.02143354791973310+01/ 

26630 

DATA  X<  3)  / 7.44160184504509300+00/ 

26640 

DATA  X<  4)  / 5.30709430617819270+00/ 

26650 

DATA  X<  5)  / 3.63401350291324620+00/ 

26660 

DATA  X<  6)  / 2.33106523030524500+00/ 

26670 

DATA  X<  7)  / 1.34479708246092630+00/ 

26680 

DATA  X<  3)  / 6.41833583695672960-01/ 

26690 

DATA  X<  9)  / 2. 01003459981210460-01/ 

26700 

DATA  X(1'0)  / 8.05943591720528330-03/ 

26710 

DATA  XSCH  1)  /0. 1 98333 1724 8562 170D  03/ 

26720 

DATA  XSCH  2)  /O. 104343885353116500  03/ 

26730 

DATA  XSQ ( 3)  /O. 553774380201781700  02/ 

.26740 

DATA  XBCH  4)  /O. 231652499746689900  02/ 

26750 

DATA  XSCH  5)  /O . 13206054139355800D  02/ 

26760 

DATA  XSCH  6)  /O ,543386510793304400  01/ 

26770 

DATA  XSCH  7)  /O. 100347919299542000  01/ 

' 26780 

DATA  XSCH  B>  /O . 412020953873336*00  00/ 

26790 

DATA  XSCH  9)  /O . 404023909244180700-01/ 

26800 

DATA  XSCH  10)  /0 . 649545073035383900-04/ 

26310 

C POSITIONS  AND  WEIGHTS  FOR  4-TERM  SUM  FOR  AIRY  FUNCTIONS 
268*^0 

DATA  U(ll)  / 4.7763903057577263D-05/ 

26830 

DATA  W( 12)  / 4.9914306432910959D-03/ 

26340 

DATA  U ( 13 ) / 8.61698469933403120-02/ 

26850 

DATA  UK  14 ) / 9.0079095845981 102D-01/ 

26860 

DATA  X<11>  / 3. 9 19832955 4 45509 ID +00/ 

26370 
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00670 

DATA  X(  12)  / 1 .6915619004823504D4  00/ 

r 

26880 

9j 

00680 

DATA  X ( 13)  / 5. 02755324672630 18D -01/ 

i 

26890 

• . 

00690 

DATA  X( 14 ) / 1 #92470605620156920-02/ 

© 

26900 

00700 

DATA  XSQ<11)  /0.15365090393596670D  02/ 

^ . 

26910 

• 

00710 

DATA  XSQ< 12)  /0.20613016631634610D  01/ 

26920 

c» 

00720 

DATA  XSCK13)  /0.2527629164866B180D  00/ 

• 

26930 

00730 

DATA  XSCK14)  /0 . 370449340277899800-03/ 

♦ 

26940 

00740 

C 

POSITIONS  AND  WEIGHTS  FOR  2-TERM  GUM  FOR  AIRY  FUNCTIONS 

26950 

o- 

00750 

DATA  WC15)  / 9.60072805957735040-01/ 

/ 

26960 

00760 

DATA  U ( 16 ) / 3. 19271940422639530-02/ 

c* 

26970 

• 

00770 

DATA  X< 15)  / 3.68006018661530440-02/ 

26980 

00780 

DATA  X( 16 ) / 1 .05924693821123730100/ 

26990 

00790 

DATA  XSQ( 15)  /0. 135423429771110700-02/ 

o 

27000 

00800 

DATA’  XSGK16)  /0. 112200407610908100  01/ 

27010 

♦ 

00810 

IF ( DX . LT . -5 . 0D0 ) GO  TO  100 

27020 

00820 

NEEDBI= . FA^SE • 

♦ 

27030 

00830 

IF<DX.GT .3.700)  GO  TO  200 

27040 

«* 

00840 

C 

THIS  ROUTE  FOR  SMALLX,  USING  POWER  SERIES.  ‘ 

27050  ' 

00950 

C 

INITIALIZE 

«> 

27060 

00860' 

10 

XS  = DX*DX 

27070 

00870 

• XCUBE  = XS  *DX 

• 

27080 

00880 

XS  = XS  *0.5D0 

• 

27090 

00890 

DF  = Cl 

27100 

• 

00900 

DFP  = C1#XS 

27110 

00910 

DG  = C2*DX 

* 

27120 

00920 

DGP  = C2 

27130 

♦ 

00930 

A I SUM  = DF  - DG 

27140 

00940 

AIPSUM  = DFP  - DGP 

• 

27150 

00950 

di  = nr  + dg 

27160 

• 

00960 

DIP  = DFP  + DGP 

27170 

00970 

FJM2=-2.0D0 

27180 

00980 

20 

FJM2 =FJN2 1 3. 0D0 ‘ 

27190 

rjMi-rjrirMOiJE 

27200 


oovvo 


-47- 


1HIS  PAGE  IS  BEST  QUAIrtfY  PRACtlCAALi  ^ 


\ 

¥ 

A 

••i 

■\ 


01000 

FJ-FJM1+0NE 

.e 

27210 

01010 

FJP1=FJ+0NE 

27220 

t 

01020 

FJP2=FJP1+0NE 

27230 

01030 

RATIO  = XCUBE/FJ 

t 

27240 

01040 

DF  = DF*RATI0/FJM1 

27250 

c 

01050 

DFP  = DFP*RATI0/FJP2 

27260 

01060 

DG  = BG*RATI0/FJP1 

c 

27270 

01070 

IiGP  = DGF*RATI0/FJM2 

27280 

c 

01080 

BI  = BI  + (BF+BG) 

27290 

01090 

BIP  = DIP  f (DFP+BGP) 

c 

27300 

01100 

IF  < NEEBBI ) GO  TO  80  ■ 

27310 

c 

OHIO 

AISUM  = AISUM  + (BF-BG) 

27320 

01120 

AIPSUM  = AIPSUM  + (BFP-BGP) 

L 

27330 

01130 

C 

CONVERGENCE  TEST 

27340 

C 

01140 

80 

IF < BABSIDF ) .GT . 1 .0D-16)  GO  TO  20 

27350 

01150 

C 

CONVERGENCE.  COMPUTE  FUNCTIONS 

c 

27360 

01160 

99 

BI  = R00T3*BI 

27370 

c 

01170 

BIP  = R00T3*BIP 

27380 

01180 

C 

THIS  RETURNS  IF  X(IS  BETWEEN  3.7  AND  8.0*  i 

c 

MORE 

27390 

01190. 

C 

ACCURATE  VALUES  OF  AI  ANB  AIP  HAVE  ALREADY 

SIAN 

27400 

c 

01*200 

C 

INTEGRATION 

• 

.27410 

01210 

/ 

IF (NEEBBI ) RETURN 

1 

• 

27420 

01220 

AI  = AISUM 

27430 

01230 

AIP  = AIPSUM 

27440 

01240 

RETURN 

27450 

01250 

c 

GAUSSIAN  INTEGRATION  FOR  LARGE  NEGATIVE  X 

27460 

'C 

01260 

100  DROOTX  = DSORT<-DX> 

27470 

01270 

R00T4X  = DSQRT (DROOTX) 

c 

27480 

01280 

DZCTA  = ~.6666666666666667*DX*DR00T:< 

27490 

c 

01290 

BARG  = DZETA  - . 7S33931633974433 

27500 

01300 

SUMR  = DZERO 

c 

27510 

01310 

SUMI  * DZERO 

27520 

c 

01320 

SUMRP  = DZERO 

C 

o 

9 

• 

© 

0 
© 
© 

y 

*> 

y 

1 

* 

* 

% 

* 


1 1 
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r 

i 


oiiio' 

SUM IP  = DZERO 

£ 

27540 

2 

01340 

C TEST  TO  SEE  HOU  MANY  TERMS  ARE  NEE 

1 

r 

01350 

27550 

. IF<HX.LT. <-2C0.D0) ) GO  TO  140 

01360 

27560 

IF<DX.LT.  ( -15.H0)  ) GO  TO  130 

r 1 

27570 

01370 

C THIS  CASE  FOR  DX  BETWEEN  -5.0  AND 

27580 

r 

01380 

LIML0=1 

01390 

27590 

LTMHI=10 

27600 

01400 

GO  TO  149  ' 

01410 

27610 

C THIS  CASE  FOR  DX  BETWEEN  -15.0  AND 

01420 

27620 

130  LIML0=11 

• 

27630 

01430 

LIMHI=14 

• 

01440 

27640 

GO  TO  149 

01450 

27650 

C THIS  CASE  FOR  BX.LT.-200. 

27660 

01460 

140  LIML0=15 

♦ 

01470 

27670 

LIMHI=16 

01480 

27680 

149  ZETASCl=DZETA**2 

' 

27690 

01490 

DO  150  K=LIMLQ « LIMHI 

- 

01500 

27700 

TERMR=U  < K > / < < ZET ASQ+XSG (K)  ) > 

01510 

27710  v 

SUMR  = SUMR  + TERMR 

27720 

01520 

TERMR=TERMR*X<K> 

01530 

27730 

3UMI*SUMI+TERMR 

27740 

01540 

TERMR=TERMR*X<K> 

27750 

01550 

SUMRP=3UME'r'+TERMR 

<r 

01560 

27760 

150  SUMIP=SUMIP+TERMR*X(K> 

27770 

01070 

SUMR*  < SUMR:»ZETASa-fSUMRP ) tZETASQ 

27780 

01580 

TEMP*SUMI-rZETASQ 

<r- 

01590 

27790 

SUMI * C TEMP-*  SUMIP  > *DZETA 

01600 

27800 

SUMRP*SUMRP*DZETA 

27810 

01610 

SUM  I(: -=SUM  IP-TEMP 

i+ 

01620 

27320 

C FORM  AIRY  FUNCTIONS 

01630 

27030 

196  S = DSIN(DARG) 

27040 

01640 

CO  = DCOG(OARG) 

- 

.01650 

27050 

RATIO  = RTPI2/R00T4X 

27060 
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U1660 

01670 


01680 


01690 


01700 


01710 


01720 


01730 


01740 


01750 


01760 


01770 


01780 


01790 


01800 


01810 


01320 


01330 


01840 


01850 


01860 


01870 

01880 


01890 


01900 


01910 


01920 


01930 


01940 


01950 


01960 


01970 


01980 


AI  = RATIO* ( COTSUNR  f S*SUMI>  . 

27870 

81  - RATI0*(C0*5UHI  - S*SUMR> 

27880 

SUMRF’=SUMRP  FSUMRP 
27890 

RATIO  - -.25D0/DX 
.27900 

FACTOR  = -RTPI2*R00T4X 
27910 

AIP  = F:ATIO*AI  - DROOTX*DI  + FACTOR* < CO*SUMRP+S*SUMIP ) 
27920 

HIP  = RATIO*BI  + DROOTX*AI  t-  FACTOR* (CO*SUMIP-S#SUMRP > 
27930 

RETURN 

27740 

C GAUSSIAN  INTEGRATION  FOR  LARGE  POSITIVE  X 
27950 

200  UROOTX  * DSQRT(DX) 

27960 

DZETA  = ,6666666666666667*DX*DR00TX 
27970 

EFAC  = DEXP< -DZETA) 

27930 

ROOT  4X  = OSORT(DROOTX) 

27990 

A I = DZERO 
28000 

BI  = DZERO 
28010 

AIP  = DZERO 
28020 

DIP  = DZERO 
28030 

IF ( DX . LT . 8 . ODO ) NEEDBI* . TRUE . 

28040 

C TEST  TO  SEE  HOW  MANY  TERMS  ARE  NEEDED  IN  GAUSSIAN  INTEGRATION 
28050 

IF(DX.GT. 15. ODO)  GO  TO  230 
28060 

C THIS  CASE  FOR  DX  BETWEEN  3.7  AND  15. 

28070 

LIML0=1 

28030 

LIMHIalO 

28090 

GO  TO  249 
28100 

C THIS  CASE  FOR  DX  GREATER  THAN  15. 

28110 

230  LIMLOail 
28120 

LIMMI=1 4 
20130 

249  DO  250  K=LIMLO» LIMHI 
28140 

DA=DZETA+X<K) 

23150 

TERrlA  = U < K ) /DA 
28160 

AI  a AI  + TERMA 
20170 

AIP=AIP+TERMA*X  < K ) /DA 
28180 

IF (NEEDS I ) GO  TO  250 
23190 
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01990 

HB-IiZETiV-XOO 

28200 

1 

02000 

TERMS  = U(K>/DD 

28210 

r 

02010 

. BI  = BI  + TERMS 

28220 

02020 

BIP=DIP+TERMB*X(K>/DB 

/• 

28230 

02030 

25 

0 CONTINUE 

28240 

02040 

C 

FORM  FUNCTIONS 

28250. 

02050 

FACTOR =RTPI*DZETA/R00T4X 

r 

23260 

02060 

RATIO  - 0.25D0/DX 

28270 

<• 

02070 

AI=AI*EFAC*FACTOR 

28280 

02080 

AIP=-(DROOTX+RATIO)*AI+RTF'I*ROOT4X*EFAC*AIP 

28290 

02090 

C 

THIS  IS  SATISFIED  ONLY  FOR  X BETWEEN  3.7  AND  3.0  IN  THESi 

SES 

28300 

02100 

C 

THE  BI  AND  DIP  ABOUT  TO  BE  COMPUTED  ARE  NOT  SUFFICIENTLY 

RATE. 

23310 

02110' 

C 

THUS  RETURN  TO  POWER  SERIES  FOR  BI  AND  DIP. 

20320 

02120 

IF ( NEEDDI ) GO  TO  10 

28330 

02130 

FACTOR =FACTOR+FACTOR 

28310 

02140  . 

BI-BIJFACTOR/EFAC 

28350 

02150 

BIP=(IiR00TX-RATI0>*BI~RTPI2*R00T4X*BIP/EFAC 

28360 

02160 

RETURN 

28370 

02170 

END 

28330 

J 
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ty  pl  r 

DIMENSION  Y ( 150) f IX <0/1 000 ) » IY<0/1000 > 

TYPE  100 

100  FORMAT < ' PLOT  TYPE=-1...N0  MORE  PLOTS'/ 

114X. '0. . .TRANSMISSION  LOSS'/ 

212Xf 'N» I . . . N-MODE  NUMBER » I-=0(MAGNITUDE)  OR  1 (PHASE)') 

1 TYPE  101 

101  FORMAT ('  ENTER  NUMBER  OF  MODES  ON  TAPE » PLOT  TYPE') 

ACCEPT  102 » MM 1 1 PL  t IMP 

102  FORMAT (100) 

NM-2*NM 

IF(IPL.EQ.-l)  STOP 
TYPE  300 

300  FORMAT ('  TYPE  TAPE  NUMBER ' ) 

ACCEPT  301 » NT 

301  FORMAT (3) 

IN=2*( IPL-1 )+IMP+l 
TYPE  103 

103  FORMAT ( ' ENTER  0 FOR  LINE  * 1 FOR  POINTS') 

ACCEPT  102 » ILPT 

TYPE  104 

104  FORMAT ( ' INPUT  NTOT r NMIN » NMAX » XM I N t XHAX f YMIM  f YHAX » MX » NY ' ) 
ACCEPT  102 » NTOT f NMIN » NMAX  f XMIN f XMAX f YMIN  f YMAX f NX f NY 
IF(IPL.NE.O)  30  TO  500 

YTEM=YMIN 
YMIN=-YMAX 
YMAX=-YTEM 
500  CONTINUE 

DX=9999./NX 
DY=9999./NY 
TYPE  130 

130  FORMAT ('  PLTP')  . 

BO  10  I— 0 f NX 
I Y < I ) =0 
IX(I)=DX*I 

10  TYPE  30? IX ( I ) f IY ( I ) 

30  FORMAT (21 10)  ' 

DO  11  I=OfNY 
. IX(I)=0 
IY( I )~DY*I 

11  TYPE  30  f IX ( I > f I Y ( I ) 

TYPE  70 

70  FORMAT ( ' PLTT' ) 

REWIND  NT 

DX=9999 . / ( XMAX-XMIN) 

DY=9999 ./(YMAX-YMIN) 

IF < IPL.NE .0)  GO  TO  200 

DO  12  1=1 f NMIN 

READ  (NTrlOS)  XfY<1) 

105  FORMAT (F13.5fFI5.6) 

IX(I)=DX*(X"XMIN) 

Y ( 1 ) =-Y ( 1 ) 

12  I Y ( 1 )=DY* ( Y ( 1 ) -YMIN) 

GO  TO  201 

200  DO  202  I=ltNMIN 

READ  ( NT »203 ) Xf (Y( J) f J=1 fNM) 

203  FORMAT (h 10 . 3 > 1P6E10 . 3/ < 10X r 1F6E10 . 3) ) 

IX ( I ) -DXY ( X -XMIN  > 

202  IY(I)=DY#(Y(IN)-YMIN) 

201  TYPE  130 

TYPE  30 f IX < NMIN ) fIY (NMIN ) 

TYPE  70 

IFdPL.NE.O)  GO  TO  204 
DU  13  I--NMINM  fNMAX 
READ  (NT.1CJ)  XfY(1) 

.::<<!)  n<v inn 
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THIS  PAGE  IS  BEST  QUALITY  PRACTICABLE 
FROM  COPY  FURMISifiEC  TO  DD.G  

yuT^yuV'  “ ' 

13  IY.(  I >-DY*< Y < 1 )-YMIN> 

GO  TO  205 

204  DO  206  I*NHIN+1»NMAX 

READ  ( NT r 203 ) X » ( Y ( J ) > J= 1 f NM ) 

IX<I)^DX*<X-XMIN5 
06  IY(I)=DY*(Y<IN>-YMIN> 

05  IF ( ILPT .EG . 1 ) GO  TO  222 
TYPE  60 

60  FORMAT ( ' PLTL ' ) 

GO  TO  223 

222  TYPE  130 

223  CONTINUE 

DO  14  I=NMINfNMAX 

14  TYPE  30  r IX  < Z ) t IY ( I ) 

TYPE  70 

GO  TO  1 


